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The  study  of  high  field  transport  has  been  instrumental  to  the 
theory  of  many  semiconductor  devices  based  on,  e.g.,  the  Hilsum-Ridley- 
Watkins  mechanism,  impact  ionization  phenomena,  and  recently  real  space 
transfer.  It  will  become  even  more  important  as  device  sizes  approach 
submicron  dimensions.  Current  interest  in  small  devices  concerns  not 
only  scaling  down  and  VLSI  (very  large  scale  integration)  but  also 
phenomena  such  as  size  quantization,  real  space  transfer,  and  velocity 
overshoot,  which  was  recently  termed  ballistic  transport. 

A  Monte  Carlo  simulation,  including  a  pseudopotential  band  structure, 
is  chosen  for  this  study.  It  is  shown  in  this  study  that  this  method 
can  be  applied  to  both  the  steady  state  and  the  transient  state  transport 
problems. 

The  speed  enhancement  of  injecting  electrons  over  the  AlxGa^_xAs-GaAs 
hetero-barrier  is  studied.  The  results  show  that  a  narrow  "collision  free 
window"  exists  with  respect  to  parameters  such  as  the  electric  field, 
the  injection  energy,  the  external  voltage,  and  the  semiconductor  dimen¬ 
sion.  It  is  found  that  only  emitter (source)-  and  base-like  structures 
are  eligible  for  collision-free  transport;  collectors (drains)  are  not 
because  of  unavoidable  high  voltage  drops. 

The  emission  of  hot  electrons  from  the  silicon  substrate  to  the 
silicon  dioxide  is  studied.  By  a  detailed  investigation  of  the  steady 
state  transport  phenomena  in  silicon,  three  sets  of  transport  parameters 
are  found.  A  Monte  Carlo  simulation  which  includes  two  realistic  con- 


duction  bands  and  the  spatial  variation  of  the  electric  field  is  then 
performed  to  study  the  high  energy  tail  of  the  distribution  function. 

Recently,  the  validity  of  the  semiclassical  Boltzmann  transport 
equation  has  been  reexamined  and  quantum  corrections  for  the  transport 
equation  have  been  suggested.  The  quantum  aspects  of  the  transport  problem 
are  discussed.  A  quantum  Monte  Carlo  method  is  proposed  and  various 
quantum  effects  are  examined.  The  study  indicates  that  the  most  impor¬ 
tant  quantum  correction  is  the  self -energy  effect  which  amounts  to  a  quasi¬ 
particle  Boltzmann  transport  equation.  The  intra-collisional  field  effect 
is  shown  not  to  be  important  in  a  steady  state  situation. 


iii 


ACKNOWLEDGMENTS 


The  author  wishes  to  express  his  deep  feeling  of  gratitude  to  his 
thesis  advisor,  Professor  Karl  Hess,  for  his  guidance,  encouragement, 
support  and  insight  during  the  course  of  this  work. 

The  author  also  thanks  Professors  N.  Holonyak  Jr. ,  G.  E.  Stillman, 
and  B.  G.  Streetman  for  many  helpful  discussions. 

The  author  is  grateful  to  Dr.  H.  Shichijo  for  many  detailed  dis¬ 
cussions  of  the  Monte  Carlo  method  in  the  early  stage  of  the  work. 

Special  thanks  go  to  Dr.  Y.  C.  Chang  for  his  many  valuable  dis¬ 
cussions  on  the  self-energy  problem. 

The  author  thanks  K.  Brennan  for  his  help  in  computer  programming 
and  proofreading  of  the  thesis. 

For  their  cooperation  and  assistance,  the  author  is  grateful  to 
Dr.  J.  P.  Leburton,  T.  C.  Hsieh,  and  M.  Keever. 

The  author  is  especially  grateful  to  S.  Brennecke  for  her  typing  of 
the  equations  and  putting  together  the  final  version. 

The  support  provided  by  the  technical  staff  at  CSL  is  gratefully 
acknowledged.  The  use  of  the  computer  facilities  of  MRL  is  also  acknow¬ 
ledged. 

The  author  is  especially  thankful  to  his  parents  for  their  continual 
support  and  encouragement . 

Finally,  the  author  wishes  to  thank  his  newlywed  wife,  Heidi.  With¬ 
out  her  love,  patience  and  encouragement,  this  work  would  not  have  been 


possible. 


iv 

TABLE  OF  CONTENTS 

CHAPTER  Page 

1  INTRODUCTION  .  1 

2  MONTE  CARLO  SIMULATION  .  4 

2.1  Introduction .  4 

2.2  Basic  Monte  Carlo  Scheme .  5 

2.3  Steady  State  Monte  Carlo .  9 

2.4  Transient  Ensemble  Monte  Carlo .  13 

3  TRANSIENT  TRANSPORT  IN  GaAs  FOLLOWING  HIGH  ENERGY 

INJECTION .  18 

3.1  Introduction . 18 

3.2  Physical  Model  and  Method  of  Computation .  19 

3.3  Results  and  Discussions .  23 

3.4  Conclusions .  35 

4  STEADY  STATE  HIGH  FIELD  TRANSPORT  IN  Si .  40 

4.1  Introduction .  40 

4.2  Theoretical  Model  .  41 

4.2.1  Band  Structure .  41 

4.2.2  Scattering  Rate .  50 

4.2.3  Temperature  Effects . 56 

4.2.4  Collision  Broadening  Effect .  57 

4.3  Keldysh's  Formalism:  P  *  100 .  57 

4.4  Comparison  and  Discussion  of  High  Field 

Transport  Parameters .  63 


V 


CHAPTER  Page 

5  STUDY  OF  ELECTRON  EMISSION  INTO  SILICON  DIOXIDE.  ...  74 

5.1  Introduction .  74 

5.2  Summary  of  Ning's  Experiments  .  75 

5.3  Method  of  Approach .  83 

5.4  Results  and  Discussions .  89 

5.5  Conclusion .  100 

6  QUANTUM  TRANSPORT  IN  SEMICONDUCTORS .  102 

6.1  Introduction .  102 

6.2  Assumptions  of  the  Boltzmann  Equation  .  103 

6.3  Quantum  Transport  and  Quasi  Particles  .  107 

6.3.1  Introduction .  107 

6.3.2  Summary  of  Barker's  Results .  108 

6.3.3  Quasi  Particle .  Ill 

6.3.4  Quantum  Monte  Carlo .  119 

6.4  Results  and  Discussion .  121 

6.5  Summary  . .  132 

APPENDIX  1  DETERMINATION  OF  CARRIER  FREE  FLIGHT  TIME 

IN  A  MONTE  CARLO  SIMULATION  .  135 

APPENDIX  2  SUMMARY  OF  SCATTERING  MECHANISMS  IN  SILICON  ....  141 

APPENDIX  3  MATERIAL  PARAMETERS  FOR  SILICON  .  153 

REFERENCES .  155 

VITA .  163 


1 


CHAPTER  1 

INTRODUCTION 

The  study  of  high  field  transport  has  been  instrumental  to  the 
theory  of  many  semiconductor  devices  based  on,  e.g. ,  the  Hilsum-Ridley- 
Watkins  mechanism,  impact  ionization  phenomena,  and  recently  real  space 
transfer  [1,2].  It  will  become  even  more  important  as  device  sizes 
approach  submicron  dimensions.  For  example,  electron  devices  with 
lengths  of  0.1-1  pm  and  supply  voltages  of  1  volt  will  operate  at  elec¬ 
tric  fields  of  10-100  kV/cm  and  corresponding  equivalent  electron  temper¬ 
atures  [3]  of  about  5000  K.  Besides  drift,  hot  electron  diffusion  and 
related  noise  phenomena  will  be  important  under  such  conditions.  Cur¬ 
rent  interest  in  small  devices  concerns  not  only  scaling  down  and  VLSI 
(very  large  scale  integration)  but  also  phenomena  such  as  size  quantiza¬ 
tion,  real  space  transfer,  and  velocity  overshoot,  which  was  recently 
termed  ballistic  transport. 

This  thesis  addresses  real  space  transfer,  the  emission  of  hot 
electrons  over  heterojunction  barriers  as  defined  by  Hess  [3],  in  a 
generalized  sense.  In  the  past,  real  space  transfer  was  mainly  considered 
in  connection  with  the  proposed  real  space  transfer  oscillator  [1].  Here 
it  will  be  shown  that  the  effect  of  hot  electrons  crossing  heterointer¬ 
faces  has  many  applications  in  the  theory  of  semiconductor  devices.  The 
following  effects  are  related  to  the  above  concept: 


(1)  the  emission  of  hot  electrons  from  silicon  into  silicon  dioxide 

(2)  the  enhancement  of  impact  ionization  for  hot  electrons  crossing 
heterointerfaces  [4] 

(3)  injection  of  electrons  over  a  hetero-barrier  to  achieve  high 
speed  enhancements  and 

(4)  effects  such  as  real  space  transfer  noise  and  hot  electron  emis¬ 
sion  out  of  buried  channels  [5]. 

The  emission  of  hot  electrons  from  the  silicon  substrate  into  the 
silicon  dioxide  [6]  has  been  known  for  some  while  and  is  based  on  the 
same  physical  concept  as  real  space  transfer.  Because  of  device  re¬ 
liability  questions,  this  latter  problem  has  received  much  attention  in 
the  past.  The  emission  of  electrons  into  SiC^  is,  of  course,  detrimental 
to  device  performance.  Thus,  real  space  transfer  can  be  used  to  con¬ 
struct  devices  but  also  appears  as  a  "nuisance"  in  device  operation. 

The  same  is,  of  course,  true  for  impact  ionization. 

Impact  ionization  phenomena  are,  on  the  one  hand,  a  limiting  factor 
in  the  reduction  of  the  device  size,  but  on  the  other  hand  represent  an 
essential  mechanism  in  the  operation  of  semiconductor  devices  such  as 
avalanche  photodiodes.  Of  special  interest  to  us  is  the  ratio  of  the 
electron  (a)  and  hole  (3)  ionization  rates.  To  minimize  the  noise  [7,8], 
a  large  a/3  (or  3/a)  ratio  is  highly  desirable.  Unfortunately,  most  III-V 
semiconductors,  including  alloys  for  long  wavelength  detectors,  exhibit 
values  of  a  nearly  equal  to  that  of  3.  In  order  to  create  an  artificially 
large  a/3  (or  8/a)  ratio,  use  can  be  made  of  the  difference  in  the  band- 
edge  discontinuities  for  electrons  and  holes  in  a  multilayered 


heterostructure.  This  is  really  a  concept  of  high  energy  injection  of 
the  carriers  inherent  to  a  heteroj unction  structure.  This  enhancement 
effect  was  first  predicted  by  theoretical  considerations  [4]  and  then 
verified  experimentally  [9]. 

This  large  variety  of  effects  could  not  be  covered  within  one  thesis 
if  they  could  not  all  be  understood  and  quantitatively  calculated  by  the 
same  approach:  A  Monte  Carlo  simulation  inclu  g  a  pseudopotential 
band  structure.  The  Monte  Carlo  simulation,  b  d  on  the  semiclassical 
Boltzmann  transport  equation,  has  been  chosen  .  ->*  .his  study.  In  Chapter 
2,  we  discuss  the  Monte  Carlo  methods  for  the  steady  state  and  the  tran¬ 
sient  state  simulation.  In  Chapter  3,  we  report  the  results  of  transport 
studies  of  GaAs  and  related  heterostructures.  In  Chapter  4,  the  band- 
structure  dependent  high  field  transport  results  for  silicon  are  discussed. 

In  Chapter  5,  the  emission  of  electrons  into  the  SiO^  from  a  silicon 
substrate  is  studied.  In  this  study,  the  high  energy  tail  of  the  hot 
electron  distribution  function  is  examined.  A  better  set  of  high  field 
transport  parameters  is  obtained  in  this  study,  which  is  a  big  step  for¬ 
ward  toward  a  better  understanding  of  high  field  transport  in  semi¬ 
conductors. 

Recently,  the  validity  of  the  semiclassical  Boltzmann  transport  equa¬ 
tion  has  been  reexamined  and  quantum  corrections  for  the  transport  equa¬ 
tion  have  been  suggested  [10,11].  In  Chapter  6,  we  discuss  the  quantum 
aspects  of  the  transport  problem.  A  quantum  Monte  Carlo  method  is  proposed 
and  various  quantum  effects  are  examined. 
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CHAPTER  2 

MONTE  CARLO  SIMULATION 


2. 1  Introduction 

The  theoretical  basis  of  high  field  transport  studies  has  been 
mainly  the  semiclassical  Boltzmann  transport  equation.  Depending  on 
the  nature  of  the  transport  problem,  different  numerical  approaches 
have  been  developed  in  solving  this  complicated  nonlinear  integro- 
differential  equation.  Among  the  many  approaches,  a  Monte  Carlo  simu¬ 
lation  has  been  the  most  successful.  It  easily  includes  all  the  scat¬ 
tering  mechanisms,  allows  the  inclusion  of  more  than  one  realistic 
energy  band  [12,13],  can  study  both  transient  and  steady  state  phenomena, 
and  has  the  power  to  extend  beyond  the  basis  of  the  classical  Boltzmann 
transport  equation  to  include  quantum  corrections  such  as  self-energy 
corrections,  collision  broadening  and  the  intracollisional  field  effect. 
Due  to  these  advantages,  the  Monte  Carlo  method  was  chosen  for  this 
thesis  study.  This  chapter  is  devoted  to  the  discussion  of  the  Monte 
Carlo  method.  As  for  the  validity  of  the  semiclassical  transport  equa¬ 
tion,  detailed  discussions  can  be  found  in  Chapter  6. 

The  Monte  Carlo  method  is  a  method  of  computer  simulation  of  a  phy¬ 
sical  system.  It  gets  its  name  from  the  use  of  random  numbers  to  simu¬ 
late  statistical  processes  in  order  to  numerically  generate  the  underlying 
probability  distribution  of  a  system.  The  application  of  simulation 
techniques  in  theoretical  physics  and  engineering  serve  two  main  purposes: 
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To  bypass  the  mathematical  difficulties  in  solving  complicated  system 
governing  equations  and  to  attain  deeper  insight  into  the  microscopic 
physical  mechanisms  that  the  governing  equations  describe.  This  method 
has  found  wide  application  in  the  field  of  statistical  physics  [14]  and 
was  first  introduced  to  the  field  of  semiconductor  transport  by 
Kurosawa  [15]  in  1966.  The  pioneering  work  of  Fawcett,  Boardman  and 
Swain  in  1970  [16],  which  successfully  reproduced  the  Gunn  effect  in 
GaAs  in  their  Monte  Carlo  simulation,  was  encouraging,  although  some  de¬ 
tails  of  the  band  structure  were  still  unknown  to  them.  Many  advances 
have  been  made  in  the  work  of  Monte  Carlo  simulation  since  then,  for  ex¬ 
ample,  by  the  Italian  group  [17-20],  the  French  group  [21-23],  etc.  The 
most  recent  advance  of  the  Monte  Carlo  method  has  been  developed  by 
Schichijo  and  Hess  who  included  a  realistic  band  structure  [12], 

In  Section  2.2,  we  discuss  the  basic  scheme  of  the  Monte  Carlo  simu¬ 
lation.  In  Sections  2.3  and  2.4,  the  steady  state  Monte  Carlo  method  and 
the  transient  ensemble  Monte  Carlo  method  are  discussed  respectively. 

2. 2  Basic  Monte  Carlo  Scheme 

There  are  two  main  ingredients  in  a  Monte  Carlo  simulation,  the  band 
structure  and  the  scattering  rates.  The  standard  Monte  Carlo  method  that 
has  been  used  by  most  researchers  incorporates  in  it  an  analytic  band 
structure,  i.e.,  a  parabolic  E(k)  relation  including  at  most  a  nonpara- 
bolicity  factor  which  accounts  partly  for  the  nonparabolic  nature  of  the 
band  structure,  and  scattering  rates  which  are  calculated  by  the  Golden 
rule,  the  effective  mass  and  the  Born  approximation  [24].  It  was  Shichijo 
and  Hess  [25]  who  first  included  the  realistic  band  structure  calculated 
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by  Che  empirical  pseudopotential  method  in  the  Monte  Carlo  simulation. 
Based  on  their  pioneering  work,  I  have  refined  and  improved  the  method 
to  include  two  conduction  bands. 

As  for  the  scattering  rates,  we  follow  the  conventional  approach 
but  modify  the  effective  mass  density  of  states  to  behave  like  the  true 
density  of  states  in  high  energies.  One  important  point  to  make  is 
that  although  scattering  rates  are  calculated  by  the  conventional  method 
to  fit  to  the  experimental  results,  the  incorrect  nature  of  the  Golden 
rule,  i.e.,  the  energy  conservation  5-function,  in  the  strong  scattering 
regime,  is  improved  in  the  simulation  by  always  selecting  the  scattering 
final  states  in  a  possible  energy  range  which  smears  out  the  sharp 
6-function  structure.  We  discuss  in  Chapter  6  the  quantum  transport  in 
semiconductors  and  a  field  theoretic  approach  to  calculate  the  high 
energy  scattering  rates. 

The  basic  Monte  Carlo  scheme  is  shown  in  Figure  2.1.  An  electron 
is  started  in  a  state  ka(i)  whose  energy  is  calculated  via  the  band  struc¬ 
ture  subroutine.  Through  the  energy  dependent  scattering  rates,  the  free 
flight  time  T(i)  of  the  electron  is  determined  by  solving 

r  =  exp 

where  r  is  a  random  number  uniformly  distributed  between  0  and  1.  The 
details  of  how  the  free  flight  time  should  be  determined  are  discussed 
in  Appendix  1.  From  the  equation  of  motion  of  a  Bloch  state  under  an 


T(i) 


dt 


T(E(k(t))) 


(2.1) 


Pseudo 
Scattering 


fRandomV 

NNumber/ 


Select 
Final 
'  State 
t 

ka(i+l)  ■ 


'  Band  ' 
Structure/ 


Figure  2.1:  Flow  chart  of  a  Monte  Carlo  simulation. 
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applied  external  field,  the  final  state  of  the  free  flight  is  calcu¬ 
lated  as 

k^i+l)  =  ka(i)  +  T(i).  (2.2) 

Again  from  the  band  structure,  the  energy  of  the  state  k^(i+l)  is  ob¬ 
tained.  A  random  number  is  then  generated  to  determine  what  scattering 

is  to  take  place.  If  it  were  a  pseudo  scattering,  the  electron  starts 

~¥  -► 

with  ka(i+l)  =  k^d+l)  again  from  the  beginning  of  the  loop.  If  it  were 
not  a  pseudo  scattering,  a  final  state  according  to  the  nature  of  the 
particular  scattering  mechanism  is  determined  through  random  numbers  and, 
of  course,  the  band  structure.  Starting  with  this  selected  new  state, 
the  electron  enters  the  beginning  of  the  loop.  This  is  the  basic  pro¬ 
cedure  of  the  Monte  Carlo  simulation. 

The  essential  idea  of  the  simulation  described  above  is  to  follow 
the  trajectory  of  the  electron  according  to  the  electronic  band  structure 
and  the  scattering  mechanisms,  and  accumulate  in  its  paths  the  relevant 
information  of  physical  quantities  like  the  drift  velocity,  the  average 
energy,  the  scattering  mean  free  path,  the  scattering  mean  free  time, 
the  distribution  function,  etc.  It  has  been  proven  that  the  distribution 
function  obtained  from  the  Monte  Carlo  simulation  solves  the  Boltzmann 
equation  [16,20]  over  a  long  simulation  path.  The  criteria  which  ensure 
that  the  observables  obtained  from  the  simulation  converge  to  the  true 
values  after  a  "long”  simulation  path  (or  time)  are  essential  to  the 
Monte  Carlo  simulation  and  are  discussed  in  Section  2.3  and  Section  2.4 
for  the  steady  state  and  the  transient  state  cases  respectively. 
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2.3  Steady  State  Monte  Carlo 

In  the  steady  state  simulation,  the  semiconductor  is  assumed  to 
be  homogeneous,  infinitely  large  and  under  a  uniform  applied  field  such 
that  the  real  space  structure  of  the  material  is  of  no  concern.  The 
trajectory  of  the  electron  is  traced  in  momentum  space  only  and  a  dis- 
tribution  function  f(k)  depending  only  on  the  crystal  momenta  is  ob¬ 
tained.  For  device  lengths  greater  than  several  microns  and  a  moderate 
spatial  variation  of  the  electric  field  in  the  device,  the  steady  state 
Monte  Carlo  simulation  can  be  applied. 

In  principle,  all  observables  of  interest  can  be  calculated  once 
f (k) ,  the  distribution  function,  is  obtained.  But  because  of  numerical 
problems,  it  is  very  costly  to  calculate  the  distribution  function  to  a 
satisfactory  degree  of  accuracy  such  that  all  relevant  quantities  can  be 
obtained  through  easy  numerical  integration  or  differentiation  of  the 
distribution  function.  To  bypass  such  numerical  difficulties,  estimators 
are  devised  to  accumulate  relevant  information  of  various  relevant  ob¬ 
servables  [26,27]. 

The  field  dependence  of  the  drift  velocity  of  a  semiconductor  is  the 
most  important  information  for  the  design  of  a  semiconductor  device. 

The  velocity  estimator  is  also  one  of  the  most  important  concerns  in  a 
steady  state  Monte  Carlo  simulation.  By  properly  choosing  a  complementary 
pair  of  velocity  estimators,  the  reach  of  the  steady  state  condition  and 
the  convergence  of  the  estimators  can  be  easily  tested.  The  two  kinds 
of  velocity  estimators  described  below  are  based  respectively  on  (1)  the 
energy  gained  from  the  external  field  during  free  flights  and  (2)  the 
energy  lost  through  phonon  emission  during  the  scatterings. 


The  first  velocity  estimator  has  been  used  by  many  authors  and 


is  given  as 


v . 
J 


^  *  3k“  dkz  =  5k7  ^f"^’ 


zi 


(2.3) 


where  the  subscript  j  refers  to  the  jth  valley,  k  and  k  ^  are  the 

initial  and  final  states  in  the  field  direction,  z,  of  the  free  flight, 

Kj  is  the  total  length  of  the  k-space  trajectory  and  the  summation  is 

-► 

over  all  the  free  flights.  ,  the  k-space  trajectory,  is  simply  re¬ 
lated  to  the  total  free  flight  time  in  the  valley  by 


r_  eF  _ 

K .  ~  T .  , 
J  "fi  3 


(2.4) 


according  to  the  equation  of  motion.  Thus  Equation  2.3  simplifies  to 


v . 
J 


E(Ef-E.) 

eFT* 

J 


ZAx 

T. 

3 


(2.5) 


where  Ax  is  the  distance  traveled  in  real  space  during  the  free  flight 
and  is  related  directly  to  the  energy  gained  from  the  field  as 


Ef  -  Ej,  =  AE  =  eFAx.  (2.6) 

So  the  first  drift  velocity  estimator  is  based  on  the  energy  gained 
from  the  field  during  the  free  flight,  which  is  directly  related  to  the 
free  flight  length  in  real  space. 

The  second  velocity  estimator  is  based  on  the  idea  of  power  bal¬ 
ance  in  a  steady  state  condition,  i.e.  ,  the  power  input  has  to  balance 
the  power  output.  The  power  balance  equation  for  one  electron  is  given 


LI 


< 


9E 

3t 


> 


eF  •  v 


(2.7) 


where  the  quantity  on  the  left  hand  side  is  the  energy  dissipation 
rate  due  to  the  phonon  scattering  events.  This  energy  dissipation  rate 
can  be  directly  accumulated  in  the  Monte  Carlo  simulation  and  hence  the 
drift  velocity  can  be  calculated  as 


The  second  velocity  estimator,  which  was  developed  in  this  work, 
is  complementary  to  the  first  drift  velocity  estimator  and  has  served 
as  the  best  convergence  test  for  the  steady  state  Monte  Carlo  simulation 
for  our  purposes. 

Depending  on  the  interests  in  the  transport  problems,  estimators 
for  different  observables  can  be  set  up  according  to  their  physical 
meanings.  For  example,  the  energy  relaxation  time  can  be  obtained  as 


Te(<E>) 


<E>  -  j  kT 

<  -  |!  >~ 
3t 


(2.9) 


where  <E>  is  the  average  energy  of  the  electron  and  is  of  course  a  func¬ 
tion  of  the  electric  field.  The  carrier  temperature  can  be  defined  and 
then  obtained  as 


<E> 

h 


(2.10) 


The  word  "defined"  means  that  the  electron  temperature  might  not  imply 
a  Maxwellian  distribution. 
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One  last  important  high  field  quantity  of  interest  to  be  discussed 
is  the  impact  ionization  coefficient.  This  estimator  is  given  as 


n , 

a  =  -iSE 
d 


(2.11) 


where  a  is  the  ionization  coefficient  defined  as  the  reciprocal  of  the 
average  distance  during  which  an  impact  ionization  event  occurs,  n^^ 
is  the  total  number  of  impact  ionization  events  and  d  is  the  total  dis¬ 
tance  traveled  by  the  carrier  in  the  field  direction  during  the  simula¬ 
tion.  In  a  steady  state  condition,  d,  the  total  distance  in  the  field 
direction  can  also  be  expressed  as 


<lc>  n  , 
f  sc 


(2.12) 


where  <lf>  is  the  mean  free  path  in  the  field  direction  between  scat¬ 
terings  and  ngc  is  the  total  number  of  scatterings.  A  simple  substitu¬ 
tion  for  d  in  a  gives 


imp 


<lc>  n 
f  sc 


(2.13) 


and  hence 


imp 


<1f> 


sc 


(2.14) 


Since  <1^>  can  be  accurately  obtained  for  a  few  hundred  scatterings  in 
the  simulation,  a  better  way  to  obtain  a  is  to  plot  the  number  of  impact 
ionization  events  as  a  function  of  the  number  of  phonon  scattering  events 
and  get  from  the  slope  of  the  curve  the  impact  ionization  coefficient  a. 
The  standard  deviation  can  also  be  estimated  from  the  plot. 
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2 . 4  Transient  Ensemble  Monte  Carlo 

As  technology  advances,  device  sizes  are  getting  smaller.  FETs 

O 

with  gate  lengths  less  than  1000  A  are  being  made  and  peculiar  hot 
electron  effects  have  been  observed  [28].  The  conventional  current 
equations  and  continuity  equation  no  longer  serve  the  purpose  in  solving 
the  transport  problem  of  such  small  devices.  One  expects  that  transient 
transport  behavior  like  the  velocity  overshoot  phenomenon  [29]  emerges 
in  such  small  devices.  To  approach  the  submicron  transport  problems 
with  a  Monte  Carlo  technique,  one  sees  immediately  that  the  real  space 
structure,  namely  the  geometry  of  the  device,  is  of  importance.  It  is 
just  like  solving  differential  equations:  initial  conditions  and  boundary 
conditions  play  the  most  important  role.  What  we  like  to  discuss  in  this 
section  is  not  the  initial  or  boundary  conditions  for  the  specific  prob¬ 
lems,  but  the  general  Monte  Carlo  techniques  for  solving  the  transient 
transport  problems  once  the  initial  and  boundary  conditions  are  specified. 

In  the  uniform  spatial  field  case,  the  basic  Monte  Carlo  scheme  for 
the  transient  simulation  is  essentially  the  same  as  that  described  in 
Section  2.2  except  that  one  also  has  to  trace  the  trajectory  of  the  car¬ 
rier  in  real  space.  The  electrons  are  released  one  by  one  from  the 
source  (cathode)  and  their  trajectories  are  stored  in  histograms  respect¬ 
ively.  The  distribution  of  quantities  of  interest  in  real  space  can  be 
obtained  by  averaging  the  ensemble  histograms  of  the  carriers.  The 
mesh  sizes  in  real  space  are  of  the  order  of  100  A  depending  on  the  ob¬ 
servables  of  interest.  The  problem  with  coo  small  a  mesh  size  is  that 


the  convergences  of  some  observables  are  hard  to  achieve.  Let  us  cite 
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an  analogy  to  explain  why  this  is  so.  If  one  has  the  grades  of  a  hun¬ 
dred  students  and  one  wants  to  plot  the  distribution  of  the  grades  (say 
from  0  to  100),  a  similar  problem  occurs  if  one  chooses  too  small  an 
interval  (say  2). 

In  testing  the  convergence  of  the  estimators,  the  best  way  to  our 
knowledge  is  to  plot  the  ensemble  average  of  the  estimators  as  a  func¬ 
tion  of  the  number  of  ensemble  carriers  and  observe  the  convergence 
from  the  curve.  Depending  on  the  observables  of  interest,  the  convergence 
curves  can  be  very  different.  For  example,  it  takes  about  500  electrons 
to  get  a  reasonable  convergence  of  the  ensemble  average  of  the  velocity 
estimator,  but  it  takes  more  than  a  thousand  electrons  to  get  a  good  es¬ 
timate  of  the  satellite  valley  populations  in  GaAs  under  fast  transient 
response  conditions. 

In  the  nonuniform  spatial  field  condition,  the  transient  Monte  Carlo 
method  needs  to  be  modified.  The  general  form  of  the  equations  of  motion 
for  a  one  dimensional  real  space  is  as  follows: 


t. 

Ak  »  ^  |  F(t)dt 


(2.15) 


and 


AE  ■  e 


F (x)dx, 


(2.16) 


where  Ak  is  the  momentum  gained  in  the  free  flight  time  from  to  t^, 
and  AE  is  the  energy  gain  during  the  acceleration.  We  see  that  Equation 
2.15  is  our  big  problem  in  that  the  field  at  time  t^<  t <  t^  depends  on 
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the  position  x(t),  the  velocity  v(x(t))  and  hence  the  band  structure 
E (k)  in  a  complicated  way.  But  if  At  =  t^-tg  is  short  enough,  which 
is  the  case  in  one  special  way  of  Monte  Carlo  simulation  (as  described 
in  Appendix  1),  we  can  proceed  as  follows.  Mathematically,  Equations 
2.15  and  2.16  can  be  written  as 


F(t)dt 


=  S?mAC  =!?(X0> 


^At  Uf  f(x0)At’ 


(2.17) 


and 

x. 

f  ^  -►  -*> 

e  F (x) dx  =  e[<{) (xn)  -  <t> ~  eF  Ax  =  E(k  )  -  E(k  )  ,  (2.18) 

J  U  i  m  1  o 

0 

where  F  is  some  unknown  mean  field  strength,  the  last  equality  in  Equation 
m 

2.17  defines  At',  y(x)  is  the  potential  profile,  and  the  relation  established 

by  the  "almost  equal  to"  sign  is  under  the  assumption  of  a  small  At. 

To  understand  the  situation  better,  let  us  examine  what  the  unknowns  and 

the  knowns  are  in  Equations  2.17  and  2.18.  The  carrier  is  now  located 

— ► 

at  the  "point"  specified  by  (tg,kg,Xg),  and  is  about  to  drift  to  the 

-V 

point  (t^.k^.x^)  in  the  interval  At  =  t^-tg.  All  quantities  related  to 
the  initial  point  are  known  but  none  related  to  the  final  point.  In 
light  of  the  last  equality  in  Equation  2.17,  one  sees  that  the  result 
of  a  carrier  drifting  in  a  time  interval  At'  under  a  uniform  electric 
field  F(Xq)  is  the  same  as  the  result  of  a  carrier  drifting  in  a  time 
interval  At  under  nonuniform  fields  provided  that  At'  is  related  to  At 


as 
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At 


F(x0) 


At'  . 


(2.19) 


F  ,  which  is  not  a  quantity  directly  related  to  the  initial  point. 


must  be  known.  We  see  that  if  we  manipulate  Equations  2.17  and  2.18, 


F  can  be  obtained  as 
m 


F 

m 


AE 

eAx 


(2.20) 


where 

AE  =  E (kg  +  -J  F (xQ) At ' )  -  E(k0)  (2.21) 

and 


Ax  =  <j> 


-1 


AE  X 

T  -  $<xo) 


-  X. 


(2.22) 


The  primary  difficulty  is  that  At'  must  be  known  to  solve  for  F  . 

m 

As  described  in  Appendix  1,  At  is  chosen  to  be  about  1/10  of 
— ► 

the  scattering  time  t/E/k^)).  As  long  as  the  choice  of  At'  does  not 
result  in  a  At  not  satisfying  the  usual  requirement  (about  1/10  of 
the  scattering  time),  we  can  proceed  in  the  simulation  to  scatter 
the  carrier  if  the  generated  random  number  r  is  such  that 


r  < 


At 


x(E(k1)) 


(2.23) 


The  reader  is  referred  to  Appendix  1  for  details.  In  the  simulation 
procedure,  the  extra  step  added  to  account  for  the  nonuniform  field 
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effect  is  simply  to  "renormalize"  the  drift  time  At  according  to 
Equation  2.19  before  determining  the  occurrences  of  the  scattering 
events.  For  a  small  spatial  variation  of  the  electric  field,  i.e., 
the  gradient  of  the  field  is  small.  At  is  not  very  different  from 
At'.  But  for  a  large  spatial  variation  of  the  electric  field,  one 
has  to  go  through  this  process. 

What  has  been  described  was  a  single  carrier  ensemble  Monte  Carlo 
technique.  There  might  be  cases  that  multicarrier  simulation  is 
necessary  and  the  inclusion  of  the  electron-electron  scattering  is  of 
importance.  In  principle,  these  can  all  be  carried  out  provided  that 
enough  money  is  invested.  It  seems  that  the  only  limitations  of  the 
Monte  Carlo  technique  in  this  respect  are  just  the  memory  and  the  speed 
of  the  computer. 

In  Chapter  6,  we  will  discuss  the  methods  of  extending  the  standard 
classical  Monte  Carlo  technique  to  include  the  quantum  effects. 
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CHAPTER  3 

TRANSIENT  TRANSPORT  IN  GaAs  FOLLOWING  HIGH  ENERGY  INJECTION 

3 . 1  Introduction 

Since  the  first  calculation  of  transient  electron  transport  in 
short  channel  FETs  by  Ruch  [29],  overshoot  phenomena  of  the  electron 
velocity  on  very  short  time  and  length  scales  have  attracted  consid¬ 
erable  attention  [39-33].  Shur  [31]  and  Shur  and  Eastman  [32]  added 
features  of  space  charge  limited  conditions  to  the  velocity  overshoot. 

They  investigated  the  initial  transient  and  called  it  the  "near  ballis¬ 
tic"  regime.  The  term  "ballistic"  is  difficult  to  define  and  is  cur¬ 
rently  applied  to  a  wide  range  of  device  parameters  and  dimensions. 
Investigations  in  References  31  and  33  indicate  that  "near  ballistic" 

O 

transport  over  larger  distances  (>  1000  A),  if  achievable,  necessitates 
injection  of  electrons  at  higher  energy.  High  energy  injection  was  also 
discussed  by  Hess  [34]  for  low  temperature  transport  free  of  scattering 

—  <4 

events  over  extremely  large  distances  L  >  10  cm  (for  electron 
energies  below  0.036  eV). 

The  criteria  for  designing  devices  and  choosing  materials  such  that 
high  transient  speeds  can  be  advantageously  achieved  using  high  energy 
injection  are  the  key  issues  of  this  study.  We  approach  the  problem  using 
a  transient  ensemble  Monte  Carlo  simulation  as  discussed  in  Chapter  2, 
which  includes  the  details  of  the  band  structure  as  calculated  by  the 
empirical  pseudopotential  method.  The  electron  is  started  at  higher 
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energy  values  (not  at  the  bottom  of  the  band  as  done  by  Ruch)  to  simulate 
injection  over  a  hetero-barrier.  The  electric  field  is  chosen  so  as  to 
be  constant  over  the  whole  distance  of  the  simulation,  in  contrast  to 
the  choice  of  the  self-consistent  field  by  Shur  [31].  The  reason  for 
this  is  the  following:  Within  the  "collision  free  window  (CFW)",  the 
electron  velocity  does  not  vary  strongly  over  the  device  length  because 
the  electrons  are  already  injected  at  high  velocities  (in  contrast  to  the 
cases  considered  by  Shur).  Therefore,  the  carrier  concentration  is 
rather  constant  and  the  electric  field  induced  by  the  carriers  is  of  minor 
importance.  We  do  not  include  the  intracollisional  field  effect  [11,35] 
because  it  would  greatly  complicate  the  computations.  In  the  CFW,  it 
would  make  only  minor  contributions  since  the  electric  fields  in  this 
window  are  small. 

With  respect  to  device  applications,  the  result  of  the  calculations 

O 

can  be  summarized  as  follows:  On  a  length  scale  of  1000  A,  emitter 
(source)-  and  base-like  structures  may  show  effects  typical  for  collision 
free  transport  (if  carefully  designed);  collector  (drain)-like  structures 
will  not. 

3 . 2  Physical  Model  and  Method  of  Computation 

As  discussed  above,  we  consider  high  energy  electrons  injected  into 
GaAs  (e.g.,  from  Al^Ga^  ^As  or  5-like  electric  fields  created  by  space 
charge  layers).  The  transition  is  assumed  to  be  abrupt,  i.e.,  the  elec¬ 
tron  gains  kinetic  energy  AE  and  forward  momentum  Ak  when  transferring 
to  the  GaAs  without  any  energy  loss. 
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As  illustrated  below,  the  calculation  of  the  self-consistent  field 
resulting  from  external  voltages  and  an  electron  redistribution  in  the 
GaAs  is  a  difficult  problem  and  depends  on  many  details,  most  impor¬ 
tantly  : 

(1)  the  boundary  conditions  of  the  injection, 

(2)  the  statistics  of  electrons  and  impurities  in  the  GaAs, 

(3)  the  level  of  injection, 

(4)  the  velocity  distribution  and  the  injection  energy. 

Let  us  assume  that  our  device  is  short  in  the  x-direction  but  rather 
wide  in  the  y-direction.  We  then  can  still  define  average  quantities 
such  as  the  density  of  electrons  in  a  meaningful  way  and  use  a  continuum 
picture  as  follows:  If  the  distribution  function  is  denoted  by  f,  we 
define  the  electron  density 


n(x) 
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or  the  current  density  as 
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dk  v.  f 
k 

—00 


(3.2) 


and  use  the  current-  and  the  continuity-equation  as  usual.  The  cal¬ 
culation  of  the  self-consistent  field,  however,  does  not  proceed  as 
simply  because  the  statistics  of  electrons  and  impurities  (donors, 
acceptors)  and  their  time  dependent  motion,  play  a  role  in  the  solution 
of  the  Poisson  equation.  As  a  consequence,  we  can  only  obtain  the 


average  electric  field 
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F  (x) 
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dy  F(x,y) 
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(3.3) 


in  a  straightforward  manner.  The  field  which  enters  the  Monte  Carlo 
calculations,  however,  (in  a  nonlinear  way)  is  F(x,y).  The  use  of  F(x) 
is  only  appropriate  under  heavy  injection  conditions  (i.e.,  large  cur¬ 
rents).  The  x-dependence  of  F(x)  is  a  function  of  the  boundary  condi¬ 
tions  and  the  average  electron  velocity.  It  can  be  obtained  for  var¬ 
ious  limiting  cases  from: 


3F (x)  _  e 
3x  £E 

o 


(3.4) 


where  NQ  represents  the  fixed  charge  (doping).  Thus  for  low  injection 
and  low  doping  density: 

F (x)  =  J  *  (3.5) 

where  V  is  the  applied  voltage.  For  heavy  injection  and  p  =  const 

F(x)  «  /x  ,  (3.6) 

and  for  F  =  const 

F(x)  *  x  .  (3.7) 

An  analysis  including  the  statistics  of  the  electron  and  impurity 
distribution  (which  appears  to  be  vital)  is  extremely  difficult.  For¬ 
tunately,  the  average  electron  velocity  does  not  sensitively  depend  on 
this  electric  field  if  the  electrons  are  injected  at  high  energies  ^nd 


if  we  restrict  ourselves  to  the  CFW,  as  will  be  seen  from  the  numerical 
results.  We  therefore  assumed  in  our  computation  the  constant  electric 
field  of  Equation  3.5.  Outside  the  CFW  and  for  high  injection  this  ap¬ 
proximation  is  weak  and  gives  only  an  estimate  of  the  approximate 
average  velocity.  The  error  can  be  estimated  from  Equations  3.5  to  3.7 
and  the  field  dependent  results  as  given  in  the  next  section. 

The  Monte  Carlo  program  used  for  the  computations  is  a  revised 
version  of  that  described  in  the  thesis  work  of  Shichijo  [25].  The 
scattering  rates  used  in  our  Monte  Carlo  model  are  identical  to  those 
used  by  Littlejohn  et  al.  [36]  for  energies  below  -  1  eV.  For  higher 
energies,  we  assume  a  proportionality  to  the  density  of  states  and,  con¬ 
sequently,  a  decreasing  scattering  rate  above  1.5  eV  [25].  However,  for 
most  results,  energies  above  1  eV  are  not  important.  But  even  for  this 
lower  energy  range  (<  1  eV) ,  the  inclusion  of  a  realistic  band  structure 
is  important.  We  chose  the  band  structure  as  calculated  by  the  empirical 
pseudopotential  method  as  described  in  [25,37].  The  energy  of  the  L-minima 
was  fixed  to  0.33  eV  and  the  energy  of  the  X-minima  to  0.522  eV.  We  use 
7000  mesh  points  of  the  E(k)  relation  calculated  by  the  empirical  pseu¬ 
dopotential  method  and  interpolate  between  these  points.  The  bandstructure 
is  not  well  simulated  by  this  model  at  very  low  energies  (few  points  only) 
and  the  results  at  low  energies  and  low  fields  are  therefore  estimated 
to  be  in  error  by  20%.  In  the  electric  field  range  of  interest  here,  the 

average  electron  energy  approaches  ~  1  eV,  where  the  total  scattering 
14  -1 

rate  is  about  10  sec  .  Our  formalism,  which  is  equivalent  to  the 
semiclassical  Boltzmann  formalism,  is  valid  in  this  range  of  energies  and 
scattering  rates  even  if  very  strict  criteria  are  applied. 


To  be  specific,  we  assume  that  Che  electrons  are  injected  from 

A1  Ga,  As  into  GaAs.  In  this  case,  one  needs  to  have  a  direct  band 
x  1-x 

gap  in  the  A1  Ga,  As,  i.e.,  x  <  0.45,  with  most  of  the  electrons  re- 
x  1— X 

siding  in  the  T  valley.  The  reason  is  that  the  electrons  transmitted 

into  the  GaAs  will  most  likely  end  up  in  the  corresponding  valleys  in 

GaAs  (because  of  conservation  of  crystal  momentum)  and  the  heavy  masses 

of  the  satellite  valleys  in  GaAs  are  undesirable  for  achieving  high 

velocity.  An  indirect  A1  Ga,  As  is  thus  undesirable  to  start  with. 

Electrons  starting  at  F  in  Al^Ga^  ^_As  will  have  (with  high  probability) 

a  wave  vector  k„  =  (k  ,0,0)  in  the  GaAs  where  E„  .  (k.)  =  AE  and  AE 
0  x  GaAs  0  c  c 

is  the  band  edge  discontinuity.  This  concludes  our  model  assumptions. 

The  numerical  results  are  given  below. 

3 . 3  Results  and  Discussions 

We  present  results  for  the  average  drift  velocity,  the  energy  dis¬ 
tribution,  the  average  number  of  scattering  events,  and  the  transit  time 
of  electrons  as  a  function  of  several  parameters. 

To  demonstrate  the  significant  effect  of  the  initial  injection 
energy,  the  average  drift  velocity  versus  distance  is  plotted  in  Figures 
3.1a  and  b  for  various  injection  energies.  The  increase  of  the  injection 
energies  is  accompanied  by  a  substantial  increase  of  the  average  drift 
velocity  of  the  electrons  as  long  as  the  injection  energy  stays  below 
the  minima  of  the  X  and  L  valleys.  Figure  3.1  clearly  demonstrates  that 
there  exists  a  CFW  with  respect  to  the  injection  energy.  Too  high  of  an 
injection  energy  works  against  the  speed  of  the  device,  as  is  seen  from 
curve  h  in  Figure  3.1a  and  curve  e  in  Figure  3.1b.  At  low  injection 
energies,  on  the  other  hand,  high  velocities  are  not  achieved. 
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Figure  3.1(a):  Average  drift  velocity  versus  the  device  length  with 
injection  energy  as  a  parameter  at  T  =  300  K.  The 
electric  field  is  applied  in  the  < 100>  direction. 
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Figure  3.1(b):  Average  drift  velocity  versus  the  device  length 
with  injection  energy  as  a  parameter  at  T  =  77  K, 
The  electric  field  is  applied  in  the  *'100> 
direction. 
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Figures  3.2a  and  b  show  the  dependence  of  the  average  drift 
velocity  on  the  strength  of  the  external  electric  field,  again  for  77  K 
and  300  K,  with  a  fixed  initial  energy  En  =  0.24  eV  (i.e.,  k  =  0.065 
2'iT/a).  By  examining  Figures  3.2a  and  b,  we  notice  that  the  peak  tran¬ 
sient  velocity  increases  with  increasing  electric  field.  However,  the 
velocity  decreases  faster  at  higher  electric  fields  as  the  length  in- 

O 

creases.  For  a  device  length  of  1000  A,  operation  in  the  CFW  regime 
requires  that  the  electric  field  does  not  exceed  ~  35  kV/cm,  which  means 
that  the  applied  voltage  must  not  exceed  0.35  V.  This  is  the  reason  for 
our  statement  that  ballistic  transport  does  not  seem  to  be  feasible  for 
collector  (drain)-like  structures  but  may  be  important  for  emitter 
(source)-  and  base-like  device  regions. 

The  physical  explanation  of  these  phenomena  is  simple.  The  electrons 
start  with  high  velocities  when  injected  into  the  central  F  valley  (small 
effective  masses)  of  GaAs  and  are  accelerated  by  the  electric  field  in 
the  forward  direction.  Those  electrons  which  survive  the  intervalley 
scattering  processes  move  up  in  energy  with  little  polar  optical  scatter¬ 
ing  and  raise  the  ensemble  average  of  the  electron  drift  velocity.  Note 
that  the  importance  of  polar  optical  scattering  decreases  with  increasing 
energy.  Since  we  investigated  electron  injection  at  high  energies,  our 
results  differ  from  the  situations  investigated  previously  [34,38],  The 
electrons  which  are  scattered  to  the  L  and  X  minima  have  a  very  low 
mobility  and  do  not  contribute  much  to  the  average  velocity.  Higher  elec¬ 
tric  fields  accelerate  the  electrons  to  the  energy  of  the  X  and  L  minima 
and  therefore  reduce  the  velocity  after  some  distance.  In  Figure  3.3,  the 


Vj  (10  cm/sec) 


Distance  (&) 

Figure  3.2(a):  Average  drift  velocity  versus  the  device  length  with 

the  electric  field  as  a  parameter  at  T  -  300  K.  The 

initial  kQ  is  given  in  units  of  2n/a,  where  "a"  is  the 

lattice  constant.  The  injection  energy  corresponding 

to  k  is  E  •  0.24  eV. 
o  o 
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Figure  3.2(b):  Average  drift  velocity  versus  the  device  length  with 

the  electric  field  as  a  parameter  at  T  =  77  K.  The 

initial  kQ  given  in  units  of  2n/a,  where  "a"  is  the 

lattice  constant.  The  injection  energy  corresponding 

to  k  is  E  =  0.24  eV. 
o  o 
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percentage  of  unscattered  electrons  is  plotted  as  a  function  of  the 
transit  length  for  different  field  strengths  at  two  temperatures  with  a 
fixed  injection  energy.  If  we  compare  Figures  3.3a  and  b  with  the  cor¬ 
responding  Figures  3.2a  and  b,  we  see  that  the  high  average  drift 
velocity  is  a  direct  result  of  the  unscattered  electrons,  which  verifies 
the  above  explanations.  Figures  3.4a  and  b  show  the  average  number  of 
scattering  events  as  a  function  of  distance.  These  plots  are,  of  course, 
related  to  Figure  3.3.  However,  a  direct  comparison  is  not  possible  be¬ 
cause  electrons  scattered  to  lower  energies  are  scattered  again  with 
higher  probability  than  the  high  energy  unscattered  electrons. 

Before  discussing  the  rest  of  the  results,  we  compare  the  differences 
introduced  by  a  temperature  variation.  The  general  characteristic  is 
similar  for  the  two  temperatures.  However,  because  of  the  freeze  out  of 
optical  phonons  at  77  K,  effects  of  collision  free  transport  are  stronger 
and  higher  transient  velocities  result.  In  other  words,  the  CFW  is  some¬ 
what  wider  for  the  lower  temperature. 

One  of  the  concerns  with  respect  to  device  operation  is  the  noise 
equivalent  temperature,  which  is  intrinsically  related  to  the  carrier 
diffusion  process  [39].  We  show  in  Figure  3.5  the  energy  distribution 

O  O 

of  electrons  at  two  distances,  700  A  and  1500  A,  for  a  field  strength  of 
20  kV/cm  and  an  injection  energy  of  0.24  eV.  The  dotted  curves  are  for 
300  K  and  the  full  curves  for  77  K.  The  broadening  of  the  distribution 
function  is,  under  certain  simplifying  assumptions,  roughly  related  to  the 
noise  equivalent  temperature  [39].  Although  the  noise  equivalent  temp¬ 
erature  is  in  general  different  from  the  average  electron  energy  (divided 
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Figure  3.3(a):  The  percentage  of  unscattered  electrons  as  a  function 

of  device  length  with  the  electric  field  as  a  parameter 
at  T  -  300  K.  These  curves  correspond  to  those  of 
Figure  3.2a. 
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Figure  3.3(b):  The  percentage  of  unscattered  electrons  as  a  function 

of  device  length  with  the  electric  field  as  a  parameter 
at  T  *  77  K.  These  curves  correspond  to  those  of 
Figure  3.2b. 


Average  Number  of  Scattering  Events 


Figure  3.4(a):  The  average  number  of  scattering  events  as  a  function 

of  device  length  with  the  electric  field  as  a  parameter 
at  T  *  300  K.  The  curves  here  correspond  to  those  of 
Figure  3.2a. 
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Figure  3.4(b):  The  average  number  of  scattering  events  as  a  function 

of  device  length  with  the  electric  field  as  a  parameter 
at  T  ■  77  K.  The  curves  here  correspond  to  those  of 
Figure  3.2b. 
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Figure  3.5:  Energy  distribution  at  two  distances  for  the  two  tempera¬ 
tures  T  ■  77  K  and  T  ■  300  K.  The  electric  field  strength 
is  10  kV/cm  applied  in  the  <100>  direction.  The  dashed 
curves  are  for  T  ■  300  K  and  the  solid  curves  are  for 
T  -  77  K. 
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by  the  Boltzmann  Constant),  it  is  clear  from  Figure  3.5  that  a  high 
noise  equivalent  temperature  results  even  for  the  initial  (almost  ballis¬ 
tic)  regime,  especially  at  room  temperature.  The  sharp  edges  of  the  77  K 
energy  distributions  are  due  to  unscattered  electrons.  High  fields  in¬ 
crease  the  average  energy  of  the  electrons  and  further  broaden  the 
energy  distribution,  as  shown  in  Figure  3.6.  It  is  interesting  to  see 
that  even  for  a  field  strength  of  100  kV/cm,  some  unscattered  electrons 
exist  for  energies  below  ~ 1  eV,  which  cause  the  small  spike  at  the  high 
energy  end  of  the  distribution  function. 

Finally,  in  Figures  3.7a  and  b,  the  ultimate  speed  (shortest  transit 
time)  of  short  GaAs  structures  is  plotted.  The  results  show  that  the 
transit  times  increase  with  increasing  electric  field.  This  paradoxical 
result  is  caused  by  the  high  scattering  rates  to  the  X  and  L  minima. 

Note  also  that  even  in  the  most  ideal  case,  it  will  be  difficult  to 

-  -  _  O 

achieve  a  transit  time  of  1  psec  for  a  transit  length  of  1000  A. 

3 . 4  Cone lus ions 

We  have  demonstrated  the  existence  of  high  transient  velocities  for 

electrons  injected  at  high  energies  into  GaAs.  To  achieve  high  speed, 

electrons  must  be  injected  at  energies  close  to  the  energy  of  the  L  minima 

and  must  not  be  accelerated  j,  ..xternal  fields  above  this  energy  (the  X 

minima  are  especially  detrimental  to  speed).  This  limits  the  injection 

energy  E^,  electric  fields  F,  and  external  voltages  V  to  rather  narrow 

0 

ranges.  Assuming  a  width  of  1000  A  for  the  structure  in  question,  we  ob- 
tain  optimum  velocities  (speed)  for  Eq  ~  300  meV,  F  <  35  kV/cm  and 
V  <  0.35  volts.  An  examination  of  these  values  immediately  suggests  that 


36 


Electron  Energy  (eV) 

Figure  3.6:  Energy  distribution  at  d  =  700  X  for  T  =  77  K.  The 
electric  fields  are  all  in  the  <100>  direction  and 
the  injection  energy  is  the  same  as  in  the  previous 
graphs . 
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Figure  3.7(b):  Transit  time  of  electrons  versus  the  device  length  with 
the  electric  field  as  a  parameter  at  77  K.  These  curves 
correspond  to  those  of  Figure  3.2b. 
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only  emitter  (source)-  and  base-like  structures  are  eligible  for  ballis¬ 
tic  transport;  collectors  (drains)  are  not  because  of  unavoidable  high 
voltage  drops.  The  high  velocities  are  caused  almost  entirely  by  elec¬ 
trons  escaping  the  intervalley  scattering  processes  which  decrease  the 
velocity  in  an  almost  step-like  manner.  Structures  substantially  longer 

O 

than  1000  A  will  not  show  ballistic  transport  except  for  very  low  elec- 
tron  energies  (just  below  the  optical  phonon  energy)  [34].  It  should  be 

O 

noted,  however,  that  even  for  transit  lengths  of  1000  A  and  below,  the 
noise  equivalent  temperature  can  be  exceedingly  high.  Because  of  the  high 
voltage  drops  at  collectors  and  drain  regions,  current  devices  will  hardly 
show  pronounced  velocity  improvement. 
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CHAPTER  4 

STEADY  STATE  HIGH  FIELD  TRANSPORT  IN  Si 

4.1  Introduction 

Silicon  has  been  the  most  extensively  studied  semiconductor 
material  since  the  invention  of  the  transistor.  Especially  the  transport 
properties  of  the  material,  which  are  important  in  device  operations, 
have  been  intensively  studied  both  experimentally  and  theoretically  by 
many  authors  [40-45].  In  the  past,  theoretical  transport  studies  have 
been  limited  to  intermediate  electric  field  strengths  (<50  kV/cm) 
largely  because  of  the  breakdown  of  the  effective  mass  approximation  at 
higher  field  strengths.  In  practical  device  operations,  electric  fields 
can  be  well  above  50  kV/cra,  for  example  at  the  drain  end  of  a  short 
channel  MOSFET.  In  order  to  penetrate  further  into  the  band  and  im¬ 
prove  the  understanding  of  the  high  field  transport  properties  of  the 
material,  one  has  to  abandon  the  effective  mass  approximation  and  adopt 
a  more  realistic  band  structure.  We  follow  the  scheme  developed  by 
Shichijo  and  Hess  [12]  who  first  incorporated  into  the  Monte  Carlo  simu¬ 
lation  a  realistic  band  structure  as  calculated  by  the  empirical  pseudo¬ 
potential  method.  We  have  improved  the  original  model  and  included  two 
conduction  bands  for  the  silicon  study. 

This  chapter  is  devoted  to  the  study  of  the  steady  state  high  field 
transport  properties  of  silicon.  Because  of  the  inclusion  of  two  con¬ 
duction  bands,  we  are  able  to  investigate  very  high  field  transport 


properties  (>100  kV/cm) ,  for  example,  impact  ionization.  In  Section  4.2, 
we  discuss  the  band  structure  and  the  scattering  mechanisms  in  silicon. 

In  Section  4.3,  we  show  the  results  of  the  simulation  by  using  the  con¬ 
ventional  Keldysh  formalism  for  the  secondary  pair  generation  rate  and 
compare  them  to  the  experimental  results.  In  Section  4.4,  we  examine 
the  large  difference  in  the  secondary  pair  generation  rates  calculated 
by  Keldysh's  formalism  and  Kane's  direct  pseudopotential  calculation. 

It  is  found  that  one  can  adjust  the  high  energy  phonon  scattering  rates 
and  the  secondary  pair  generation  rates  to  both  (Keldysh  and  Kane)  models 
and  still  find  a  decent  fit  of  the  theoretical  results  to  the  experimen¬ 
tal  data. 

4 . 2  Theoretical  Model 

The  model  for  Monte  Carlo  simulation  has  two  main  ingredients:  the 
band  structure  and  the  scattering  rate.  We  described  in  the  following 
sections  the  different  features  that  have  been  included  in  our  model  and 
the  advantages  it  has  over  other  models. 

4.2.1  Band  Structure 

Figure  4.1  shows  the  band  structure  of  silicon  calculated  by  the 
empirical  pseudopotential  method  [37].  Note  that  the  first  conduction 
band  and  the  second  conduction  band  intersect  at  the  X  point,  where  the 
energy  of  the  state  is  about  0.1  eV.  For  a  moderately  high  electric 
field,  the  average  energy  of  the  electrons  easily  gets  above  0.1  eV  and 
interband  transitions  occur.  Unlike  the  case  for  GaAs ,  where  the  second 
conduction  band  lies  significantly  higher,  one  has  to  include  the  second 
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Figure  4.1:  The  band  structure  of  silicon  calculated  by  the 
empirical  pseudopotential  method  [37]. 
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band  in  the  case  for  silicon.  Note,  also,  that  the  L  minimum  is  only 
1  eV  above  the  X  minimum.  Under  very  high  electric  fields  ('•100  kV/cm) 
where  impact  ionization  phenomena  are  important,  one  can  expect  that  the 
L  valleys  play  an  important  role. 

Another  transport  process  of  interest  to  us,  which  is  discussed 
in  Chapter  5,  is  the  emission  of  electrons  into  silicon  dioxide.  The 
potential  barrier  at  the  interface  of  silicon  and  silicon  dioxide  is 
about  3.1  eV.  It  is  clearly  necessary  to  look  carefully  into  the  band 
structure  and  see  where  in  the  Brillouin  zone  such  a  process  is  possible 
before  investigating  transport  at  such  high  energies.  Figure  4.2  shows 
a  cross  section  of  the  Brillouin  zone.  The  isoenergy  lines  corresponding 
to  this  cross  section  are  shown  in  Figures  4.3a  and  b  for  the  first  and 
the  second  conduction  bands  respectively.  We  see  that  the  X  valleys  of 
the  first  conduction  band  are  more  elliptic,  while  for  the  second  con¬ 
duction  band,  they  are  more  isotropic.  Note  that  there  is  only  a  small 
region  near  the  W  point  for  the  first  conduction  band  that  is  above  3  eV 
on  this  cross  section.  Figure  4.4  shows  another  cross  section  of  the 
Brillouin  zone,  and  Figures  4.5a  and  b  illustrate  the  isoenergv  lines  cor¬ 
responding  to  that  cross  section  for  the  two  bands.  Note  that,  on  this 
cross  section,  the  energies  are  always  below  3  eV  for  the  first  conduction 
band.  Thus,  it  is  essential  that  one  includes  the  second  conduction  band 
in  order  to  look  at  high  energy  processes  like  the  emission  of  electrons 
over  the  silicon-silicon  dioxide  barrier. 

Because  of  the  symmetry  properties,  we  only  need  to  calculate 
points  inside  the  irreducible  wedge  of  the  Brillouin  zone.  By  means  of 
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Figure  4.3(a):  Equal  energy  surfaces  for  the  .  t 

conduction  band  of  silicon  corresponding 
to  the  cross  section  of  Figure  4.2. 
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5(a) :  Equal  energy  surfaces  for  the  first 

conduction  band  of  silicon  corresponding 
to  the  cross  section  of  Figure  4.4. 
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the  48-fold  symmetry  operations,  one  easily  recovers  the  whole 
Brillouin  zone. 

As  described  in  Reference  25,  152  points  inside  the  irreducible 
wedge  plus  another  97  points  around  the  wedge  (necessary  for  interpo¬ 
lations)  were  directly  calculated  by  the  empirical  pseudopotential 

■+ 

method.  Table  4.1  lists  the  calculated  E(k)  relation  for  these  249 
points . 

4.2.2  Scattering  Rate 

The  scattering  mechanisms  included  in  our  Monte  Carlo  simulation 
are  the  following:  intravalley  acoustic  scattering,  X  to  X  eouivalent 
intervalley  f-scattering ,  X  to  X  equivalent  intervalley  g-scattering , 
and  X  to  L  nonequivalent  intervalley  scattering,  which  has  never  been 
considered  in  previous  works.  By  definition,  f-scattering  is  the 
scattering  of  electrons  to  any  of  the  four  neighboring  X  valleys,  for 
example,  from  (100)  valley  to  (010)  valley,  and  g-scattaring  is  the 
scattering  of  electrons  to  the.  opposite  valley,  for  example,  from  (100) 
valley  to  (-100)  valley.  We  follow  closely  the  work  by  Canal!  et  al. 
in  calculating  the  low  energy  scattering  rates  [13].  For  X  to  X  coupling, 
we  consider  several  possible  phonon  types  (3f  and  3g)  [18,42].  The  dif¬ 
ference  between  f  and  g  scatterings  is  taken  into  account  in  the  simula¬ 
tion  process  by  properly  selecting  the  scattering  final  states.  For 
intravalley  acoustic  scattering,  the  energy  exchange  between  electrons 
and  phonons  is  taken  into  account.  The  X  to  L  coupling  constants  were 
assumed  to  be  the  same  for  all  four  possible  phonons  determined  from  the 


TubLe  4.1a:  K(k)  relation  for  the  lowest  conduct  ion  band  of  Si 
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phonon  spectrum  [46-48]  and  were  found  by  fitting  the  calculated  results 
to  the  experimental  results. 

The  effective  mass  density  of  states  with  the  nonparabolicity 
factor  included  was  used  in  calculating  the  low  energy  scattering  rates. 

It  is  not  appropriate  to  simply  extend  these  scattering  rates  to  higher 
energies  because  the  effective  mass  approximation  breaks  down  at  higher 
energies  and  the  density  of  states  starts  decreasing  at  some  critical 
point.  The  intervalley  scattering  rate,  which  is  proportional  to  the 
final  density  of  states,  should  be  modified  at  higher  electron  energies. 
Therefore  we  have  calculated  the  total  density  of  states  for  the  two 
conduction  bands  for  Si.  As  discussed  in  Section  2.2,  we  modify  the  high 
energy  scattering  rates  to  behave  like  the  density  of  states.  The  de¬ 
tails  of  the  mathematics  in  calculating  the  scattering  rates  are  summarized 
in  Appendix  2. 

For  the  interband  coupling,  we  consider  the  interband  and  intraband 
scatterings  with  the  same  type  of  phonons  and  the  same  coupling  strength. 

In  other  words,  electrons  see  the  same  scattering  rates  both  in  bands  1 
and  2  and  get  scattered  equally  likely  into  possible  final  states  in  both 
bands.  This  simplifies  drastically  the  problem  with  interband  and  intra¬ 
band  scatterings  in  the  simulation.  One  can  assume  independent  coupling 
parameters  for  interband  scatterings  and  adjust  them  to  fit  to  the  experi¬ 
ment.  But  with  interband  deformation  potentials  and  interband  overlap 
integrals  essentially  unknown,  one  currently  does  not  refine  the  model  by 
adding  in  extra  parameters. 

Impact  ionization  can  be  treated  as  an  additional  scattering 
mechanism  in  the  Monte  Carlo  simulation.  In  previous  works  [49-54],  this 


secondary  pair  generation  rate  was  always  calculated  by  Keldysh's 
empirical  formula: 

f F  -  E  , 
th 
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(4.1) 


where  E  is  the  threshold  energy  for  impact  ionization,  t  (E  )  is 
th  ph  th 

the  total  phonon  scattering  rate  at  the  ionization  threshold  energy 
^th’  an<^  P  is  a  constant.  It  has  been  shown  that  the  impact  ionization 
coefficient  is  essentially  unchanged  as  long  as  P  is  much  greater  than 
one.  This  amounts  to  saying  that  impact  ionizations  occur  almost  im¬ 
mediately  for  electrons  with  energies  above  the  threshold  energy. 

Kane  has  done  a  direct  calculation  of  the  secondary  pair  genera¬ 
tion  rate  in  silicon  using  the  pseudopotential  wave  function  and  assuming 
two  particle  interaction  only  [55].  Surprisingly,  his  calculated  rate 
is  much  lower  than  the  rate  calculated  from  Keldysh's  formula.  We  have 
performed  a  full  investigation  on  the  discrepancy.  In  Section  4.3,  we 
show  the  results  of  the  simulation  following  the  conventional  Keldysh 
formalism.  In  Section  4.4,  we  compare,  in  detail,  the  results  for  three 
cases:  1.  the  rate  calculated  from  the  conventional  Keldysh  formalism 

with  P> > 1 ;  2.  the  rate  calculated  from  Keldysh's  formalism  but  with  P < 1 ; 
3.  the  rate  from  Kane's  direct  pseudopotential  calculation. 


4.2.3  Temperature  Effects 

We  assume  that  the  band  structure  is  essentially  unaltered  bv  the 


temperature  effect,  except  that  the  band  gap  is  modified.  We  further 


assume  that  the  threshold  energy  of  impact  ionization  can  be  scaled 
according  to  the  formula 
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The  formula  for  the  temperature  dependence  of  the  band  gap  E^  can  be 
found  in  [56].  With  the  above  assumptions,  the  temperature  variation 
on  impact  ionization  coefficient  can  be  investigated  by  using  Keldysh's 
formalism. 


4.2.4  Collision  Broadening  Effect 

14 

As  the  scattering  rate  goes  up  to  10  per  second,  the  energy  ot 
the  electron  is  uncertain  within  K/t  according  to  the  uncertainty 
principle.  AE  is  about  100  meV  for  1/:  higher  than  lO^"1.  In  determining 
the  final  state  after  each  scattering,  we  consider  the  states  within  a 
fixed  range  AE  of  the  true  final  state  energy  as  possible  candidates. 

We  then  randomly  choose  one  of  those  candidates  as  the  final  state,  keep¬ 
ing  r'  e  average  energy  loss  (gain)  constant.  Part  of  the  collision 
broadening  effect  should  be  taken  into  account  automatically  by  so  doing. 
In  Chapter  6,  a  better  method  taking  into  account  the  Lorentzian  shape 
of  the  collision  broadening  is  discussed.  For  all  simulations  in  Chapters 
5  and  6,  the  collision  broadening  effect  is  only  partly  accounted  for  as 
described  above. 

4.3.  Keldysh's  Formalism  :  P  =  100 

Using  Keldysh's  formula  with  P  =  100,  we  have  calculated  drift 


velocities  in  the  < 1 1 1 >  direction  for  high  electric  fields  and  compared 
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them  to  the  recent  high  field  experimental  results  by  Smith  et  al.  [18, 
57].  As  is  shown  in  Figure  4.6,  our  calculated  high  field  results  are 
in  good  agreement  with  the  experimental  results.  In  Figure  4.7,  we  show 
the  calculated  results  of  the  electron  impact  ionization  rate  at  two 
temperatures,  100  K  and  300  K.  For  300  K  our  results  seem  to  agree 
better  with  Lee  and  Logan's  results  than  with  Overstraeten  and  De  Man's 
results.  We  did  not  find  any  orientational  dependent  ionization  rate  in 
Si.  Since  the  previous  work  by  Shichijo  et  al.  [12]  did  not  find  any 
orientational  dependence  of  the  electron  ionization  coefficient  in  GaAs 
with  the  same  kind  of  approach,  it  leads  us  to  believe  that  the  orienta¬ 
tional  effect  is  connected  with  some  complicated  effects. 

As  discussed  by  Hess  et  al.  in  a  recent  paper  [61],  the  answer  to 
the  anisotropy  may  lie  in  the  inhomogeneity  of  the  electric  field.  In 
the  transient  regime,  the  statistics  of  impurities  most  probably  is  an 
important  factor  and  it  is  conceivable  that  aggregates  of,  for  example, 
Poisson  distributed  impurities,  led  to  high  local  field  fluctuations  and 
impact  ionization  enhancements  as  described  by  Shockley  [50]  and  recently 
for  superlattices  bv  Chin  et  al.  [4].  The  corresponding  enhancement  of 
the  electron  ionization  coefficient  as  described  bv  Chin  et  al.  is 
essentially  ballistic  and  may  be  a  source  of  the  anisotropy.  As  we  men¬ 
tioned  in  Section  4.2,  the  first  two  conduction  bands  were  included  in 
the  model.  We  show  some  interesting  but  important  results  in  Figure  4.8 
which  tell  us  the  important  role  played  bv  the  second  conduction  band  in 
the  transport  of  electrons  in  Si.  We  plotted  on  the  graph  both  the 
second  band  effect  and  the  average  energy  of  electrons  as  a  function  of 
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direction.  The  dashed  line  and  the  solid  line  are  the 
experimental  reults. 
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Figure  4.7:  Calculated  electron  ionization  coefficient  i 


plotted  as  a  function  of  inverse  field  strength. 
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4.8:  The  second  band  effect  and  the  average  energy  o:  electrons 

plotted  as  a  function  of  the  electric  field.  The  iert 
scale  corresponds  to  the  upper  curve  and  the  right  scale 
corresponds  to  the  lower  curve. 
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electric  field.  From  the  pseudopotential  band  structure,  the  minima  of 
the  second  conduction  band  are  located  exactly  at  X.  They  are  about 
0.1  eV  above  the  minima  of  the  first  conduction  band.  The  electrons  can 
be  scattered  to  the  second  band  once  their  energies  get  above  0.1  eV. 

The  correlation  between  the  second  band  effect  and  the  average  electron 
energy  shows  clearly  in  Figure  4.8.  The  decrease  of  the  effect  of  the 
second  band  after  100  kV/cm  is  probably  because  of  the  slower  relative 
increase  of  the  density  of  states  of  the  second  band  compared  to  the 
first  band.  Previous  Monte  Carlo  calculations  by  the  Italian  group  [13, 

62]  produced  a  drift  velocity  consistently  higher  than  the  experimental 
values.  We  believe  the  reason  for  this  lies  in  the  neglect  of  the  second 
band,  which  enhances  the  scattering  rate  of  electrons  for  energies  above 
0.1  eV.  The  scattering  percentage  due  to  X  to  L  scattering,  as  seen  from 
our  simulation,  is  significant  only  when  the  electric  field  is  above  100 
kV/cm.  The  drift  velocity  for  electric  fields  under  100  kV/cm  is  basically 
not  influenced  by  adjusting  the  X  to  L  coupling  constant.  The  impact 
ionization  rate  is  verv  sensitive  to  this  coupling  constant  because  the 
X  to  L  scattering  is  sign;:  :  i:it  wnen  the  electric  field  is  above  200  kV/em. 
For  an  assumed  i.sotr'pi  ■  n  Id  energy  in  Si  of  1.8  eV  (60],  a  coupling 
constant  of  3  :<  10  eV  cm  i  i  'stained  hv  fitting  the  experimental  data 
at  333  kV/cm. 
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4 .  4  Comparison  and  Discussion  of  High  Field  Transport  Parameters 

As  mentioned  in  the  last  section,  one  can  adjust  the  X  to  L  cou¬ 
pling  constant  and  the  impact  ionization  scattering  rate  to  fit  the  ex¬ 
perimental  data.  In  Figure  4.9,  we  show  the  scattering  rates  as  a 
function  of  energy  for  three  different  sets  of  parameters  which  all  fit 
the  experimental  results  within  the  error  of  the  calculations.  The 
scattering  rate  below  1  eV,  which  is  important  in  determining  the  low 
field  drift  velocity  (<100  kV/cm),  stays  the  same  for  the  three  cases. 

In  other  words,  the  low  field  transport  behavior  for  the  three  sets  of 
parameters  are  the  same.  The  parameters  of  Set  1  are  those  used  in 

Q 

Section  4.3  with  DVT  =  3  x  10  eV/cm,  P  =  100,  and  E  ,  =  1.8  eV.  The 
XL  tn 

g 

parameters  of  Set  3  are  based  on  D  =  1  x  10  eV/cm  and  Kane's  impact 

AL 

ionization  scattering  rate  [55].  The  parameters  of  Set  2  lie  somewhat 

g 

in  between  the  two  sets  of  parameters  and  are  D  =  2  x  10  eV/cm, 

A  Li 

P  =  0.01  and  E  ^  =  1.1  eV  for  Keldysh's  formula.  Figure  4.10  shows  the 
impact  ionization  scattering  rate  for  the  three  cases.  Note  that  for 
P  =  100,  the  ionization  scattering  rate  shoots  up  almost  abruptly  to 
very  high  values  and  dominates  over  the  phonon  scattering  rate.  For  Sets 
2  and  3,  the  ionization  scattering  rates  are  more  moderate  and  phonon 
scatterings  remain  dominant  up  to  very  high  energies  (about  2.5  eV).  In 
Figure  4.11,  we  show  the  electron  ionization  coefficient  plotted  against 
the  inverse  field  strength  for  the  three  sets  of  parameters.  Within  the 
standard  deviation  of  the  calculations,  they  all  lie  in  the  range  of  the 
experimental  results.  It  is  surprising  that  all  three  sets  of  parameters 
fit  to  the  experimental  results  in  spite  of  the  fact  that  the  behavior  of 
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Figure  4.9:  The  total  scattering  rates  corresponding  to  the  three  sets 
of  transport  parameters  described  in  the  text  plotted  as  a 
function  of  energy.  The  upper  curve  for  each  set  at  high 
energies  corresponds  to  the  impact  ionization  scattering  rate. 
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Impact  ionization  scattering  rates  for  the  three  sets  of 
transport  parameters  plotted  as  a  function  of  energy. 
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Figure  4.11:  Electron  ionization  rate  plotted  against  the  inverse  field 
strength  for  the  three  sets  of  transport  parameters. 
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the  electrons  is  quite  different.  In  the  situation  for  Set  1  parameters, 
which  is  the  one  that  has  been  used  for  all  previous  impact  ionization 
theories,  the  electrons  that  gain  energy  above  the  ionization  threshold 
impact  ionize  almost  immediately.  The  effective  threshold  energy  is 
only  a  little  higher  than  the  threshold  energy.  For  the  other  two  sets 
with  threshold  energy  equaling  the  band  gap  of  1.1  eV,  the  electron  does 
not  impact  ionize  immediately  after  its  energy  exceeds  the  threshold. 

The  phonon  scatterings  and  impact  ionization  are  then  competing  and  the 
effective  threshold  becomes  a  strong  function  of  the  electric  field  as 
shown  in  Figure  4.12.  In  the  case  of  Kane's  ionization  rate,  the  ef¬ 
fective  threshold  can  be  as  high  as  3.3  eV  for  a  field  strength  of  400 
kV/cm. 

The  calculated  drift  velocities  corresponding  to  Figure  4.11  are 
shown  in  Figure  4.13.  Notice  that  the  lower  the  high  energy  scattering 
rate,  the  higher  the  high  field  drift  velocities.  The  reason  for  this 
trend  can  be  explained  as  follows.  Intervalley  scatterings  are  random¬ 
izing  and,  hence,  randomize  the  forward  momenta  of  the  electrons.  The 
magnitude  of  the  average  forward  momentum  in  the  field  direction  is  dir¬ 
ectly  proportional  to  the  drift  velocity.  Therefore,  the  less  the 
randomizing  scattering,  the  higher  the  drift  velocities.  However,  this 
argument  is  only  partially  true  because  the  mass  of  the  electron  is  a 
strong  function  of  energy.  But  the  subtleties  are  seen  when  reexamin¬ 
ing  Figure  4.11. 

There  is  no  obvious  trend  for  the  three  sets  of  results  shown  in 
Figure  4.11.  In  fact,  there  is  a  subtle  balance  between  the  phonon 
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Figure  4.13:  High  field  drift  velocities  for  the  three  sets  of  trans¬ 
port  parameters  plotted  against  the  electric  field. 
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scattering  process  and  the  impact  ionization  process.  In  order  to 
really  understand  the  details,  one  has  to  resort  to  the  distribution 
function. 

The  strength  of  the  total  scattering  rate  (including  the  phonon 
scattering  process  and  the  impact  ionization  process)  determines  how  much 
the  distribution  function  is  "heated  up".  As  shown  in  Figure  4.14,  the 
lower  the  total  scattering  rate,  the  higher  the  electron  temperature  as 
defined  by  Equation  2.10.  In  other  words,  there  are  more  electrons  with 
high  energies  if  the  total  scattering  rate  is  lower.  But  the  larger  num¬ 
ber  of  electrons  in  the  high  energy  tail  of  the  distribution  function  does 
not  imply  a  larger  ionization  coefficient  because  there  are  two  scattering 
processes,  i.e.,  phonon  scattering  and  impact  ionization  scattering,  compet¬ 
ing  with  each  other. 

Figure  4.15  shows  the  total  energy  dissipation  rates  and  the  impact 
ionization  dissipation  rates  for  the  three  cases.  One  might  expect  that 
the  high  ionization  rate  results  in  a  high  impact  ionization  dissipation 
rate.  It  would  seem  that  the  faster  (in  time)  the  electrons  impact  ionize, 
the  higher  the  dissipation  rate.  But  Figure  4.15  shows  just  the  opposite. 

It  is  simply  because  the  total  rate  of  an  event  happening  is  given  by  the 
product  of  the  number  of  particles  available  for  the  event  and  the  in¬ 
dividual  rate  for  the  given  event.  From  Figure  4.14,  we  see  that  Set  3 
gives  the  highest  electron  temperature,  which  implies  a  larger  number  of 
electrons  in  the  high  energy  tail,  and  hence  the  highest  dissipation  rate 
as  shown  in  Figure  4.14,  although  its  total  scattering  rate  is  the  lowest. 

How  can  we  distinguish  among  the  three  sets  of  high  field  trans¬ 
port  parameters,  which  all  agree  well  enough  with  the  experiments?  As 
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Figure  4.15:  The  total  energy  dissipation  rate  and  the  impact  ionization 
dissipation  rate  plotted  against  the  electric  field  for  the 
three  sets  of  transport  parameters. 
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shown  in  Figure  4.14,  the  electron  temperature,  which  is  a  measure  of 
the  energy  distribution  (although  it  does  not  necessarilv  imply  a 
Maxwellian  form),  is  very  different  for  the  three  cases.  Therefore,  if 
one  can  measure  the  distribution  function  at  high  energies  by  some  ex¬ 
perimental  means,  one  should  be  able  to  determine  the  right  set  of 
parameters . 

In  previous  theories  of  impact  ionization  [49,  54],  the  ionization 
cross  sections  were  all  assumed  to  be  abruptly  rising  above  the  threshold 
energy.  As  clearly  shown  in  Figure  4.9,  this  is  not  necessarily  the  case. 
We  describe  in  the  next  chapter  how  Ming's  experiment  measures  the  high 
energy  tail  of  the  distribution  function  and  how  we  simulate  the  experi¬ 
ment  by  a  Monte  Carlo  method.  It  turns  out  that  Set  2  of  the  parameters 
describes  the  high  energy  tail  of  the  distribution  function  correctly  and 
the  previous  assessments  of  impact  ionization  increasing  rapidly  above  a 


threshold  of  1.5E 

g 


[56]  seems  to  be  in  error. 


CHAPTER  5 


STUDY  OF  ELECTRON  EMISSION  INTO  SILICON  DIOXIDE 

5 . 1  Introduction 

The  Si/Si07  system  has  been  studied  extensively  because  of  its 
important  role  in  silicon  device  technology.  One  of  the  most  intriguing 
hot  carrier  transport  phenomena  in  this  structure  has  been  the  emission 
of  hot  electrons  or  holes  into  the  Si07  layer.  This  phenomenon  is,  on 
the  one  hand,  an  unwelcome  effect  causing  threshold  shifts,  breakdown 
walkout  and  hence  device  instability  [63-65);  but  on  the  other  hand,  it 
has  found  application  in  electrically  programmable  memory  devices  [66]. 

In  order  to  understand  the  emission  process,  various  experiments  have  been 
constructed  to  study  this  problem  [67-73],  Among  the  many  attempts, 

Ning's  experimental  setup  has  been  the  most  interesting  one.  He  was  able 
to  measure  the  absolute  emission  probability  of  optically  induced  in¬ 
jection  electrons  coming  from  the  substrate  toward  the  Si-SiO,  interface. 
In  another  sense,  he  was  able  to  probe  the  high  energy  tail  of  the  dis¬ 
tribution  function,  which  is  extremeiv  important  in  determining  the  high 
field  transport  parameters. 

Past  models  explaining  the  emission  process  invoLved  simplifying 
assumptions,  for  example,  an  energv  independent  scattering  mean  free  path 
and  a  parabolic  band  structure.  As  discussed  in  Section  -*.2,  the  hand 
structure  is  highly  nonparabolic  at  high  energies.  Such  high  energy 
transport  has  to  consider  different  valley  tvpes  (X.L)  and  more  than  one 
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conduction  band.  Moreover,  impact  ionization  definitely  plays  an  im¬ 
portant  role.  Since  Ming's  experiments  involved  low  current  injection, 
they  are  best  suited  for  a  band  structure  dependent  Monte  Carlo  study. 

This  chapter  is  devoted  to  the  theoretical  explaining  of  Ming's 
experiment  with  two  purposes  in  mind:  To  understand  the  physics  behind 
the  emission  process  and  to  determine  the  right  set  of  high  field  trans¬ 
port  parameters  to  probe  farther  up  into  the  band.  In  Section  5.2,  a 
summary  of  Ming's  experimental  procedure  is  given.  In  Section  5.3,  we 
discuss  the  Monte  Carlo  method  for  the  simulation.  In  Section  5.4,  the 
results  and  discussion  are  given.  In  Section  5.5,  we  present  our  con¬ 


clusion. 


5 . 2  Summary  of  Ming's  Experiments 

The  devices  used  in  Ming's  experiment  were  n-channel  polysilicon- 
SiO,-Si  field-effect  transistors.  Use  has  been  made  of  the  optically  in¬ 
duced  hot-electron  injection  [69]  as  illustrated  in  Figures  5.1a  and  b.  In 
this  method,  source  and  drain  were  grounded,  a  negative  bias  was  applied 
to  the  substrate,  and  a  positive  bias  was  applied  to  the  gate.  The  advan¬ 
tage  of  this  setup  is  that  the  gate  voltage  and  the  substrate  voltage  can 
be  independently  adjusted  and  are  not  affected  by  each  other  since  the 
interface  is  pinned  to  the  ground  potential.  Electron-hold  pairs  were 
generated  in  the  substrate  by  incident  photons.  Electrons  which  diffused 
into  the  depletion  region  are  accelerated  toward  the  Si-SiO ,  Interface. 

The  majority  of  carriers  thtt  do  not  >vercome  the  interface  barrier  are 
collected  as  the  source  and  the  drain  currents.  The  carriers  that  overcome 
the  barrier  are  collected  as  the  gate  current  assuming  that  they  have  not 
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been  trapped  in  the  silicon  dioxide  layer.  By  measuring  the  gate,  source 
and  drain  currents,  the  absolute  emission  probability  could  be  determined 


by 


P 


'total 


(5.1) 


where  I  .  was  the  total  current  coming  from  the  substrate  and  was  ap- 
total  °  r 

proximately  twice  the  drain  current,  since  the  gate  current  is  orders  of 
magnitude  smaller  than  the  drain  current,  and  the  source  and  the  drain 
currents  were  equal  by  the  symmetry  of  the  device. 

The  doping  profile  in  the  silicon  substrate  can  be  represented  by 
a  Gaussian  of  the  form 

2 

N  (x)  =  N  +  C  exp(--~)  ,  (5.2) 

A  B  0  2cf 

where  x  is  the  distance  from  the  Si-SiO,  interface,  N_  is  the  background 
doping  concentration,  is  the  surface  concentration  due  to  boron  im¬ 
plantation,  and  C  is  the  characteristic  of  the  penetration  of  C^.  The 
solution  to  the  Poisson  equation  for  this  doping  profile  is  derived  in 
detail  in  the  appendix  of  Reference  73. 

Figures  5.2a  and  b  show  the  potential  profiles  and  the  electric 
field  profiles  in  the  depletion  region  for  two  substrate  voltages  with 
Ng  =  7.5  x  1015  cm  \  =  1.13  x  10*'**  cm  ^ ,  and  a  =  3.36  x  10  D  cm  (device 

15-2-9  in  Ming's  experiment).  Not  ice  that  the  potential  profiles  are 
strong  functions  of  the  spatial  coordinates  and  hence  the  electric  fields 
are  highly  nonuniform.  This  has  made  it  too  difficult  to  approach  the 
problem  by  solving  the  Boltzmann  equation  directly  [68,  75]. 


Potential  Profile  (V) 
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Figure  5.2(a):  Potential  profiles  in  the  depletion  region 


tor  device  15-2-9. 


Electric  Field  (kV/cm) 
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Figure  5.2(b):  Electric  field  profiles  in  the  depletion 
region  for  device  15-2-9. 
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The  luck_,  electron  model  proposed  by  Ming  et  al.  to  explain  the 
emission  process  is  depicted  in  Figure  5.3,  The  concept  of  a  "lucky 
electron"  was  conceived  by  Shockley  in  treating  the  impact  ionization 
problem  [50],  Electrons  located  a  distance  d  away  from  the  interface 
having  a  potential  energy  equal  to  the  emission  barrier  will  overcome  the 
barrier  and  go  into  the  silicon  dioxide  layer  if  they  encounter  no 
scatterings  at  all.  For  an  assumed  constant  scattering  mean  free  path  a, 
the  probability  for  lucky  electrons  can  be  expressed  as 

P  =  A  exp  (-  -r),  (5.3) 

where  A  is  a  normalization  constant.  This  expression  can  be  derived 
Che  same  wav  as  we  derived  the  time  of  free  flight  in  the  Monte  Carlo 
simulation  (Appendix  1). 

The  barrier  lowering  effect,  which  is  also  shown  in  Figure  5.3, 
is  due  to  Schottky  lowering  and  effective  tunneling  lowering.  According 
to  Ning,  the  barrier  height  can  be  written  in  the  form 

1/2  2/3 

if  =  3.1  eV  -  b  F  '  -  a  F  ,  (5.4) 

B  ox  ox 

where  3.1  eV  is  the  Si-SiCL  interface  barrier  [76],  F  is  the  oxide 

2  ox 

field.  The  second  terra  is  responsible  for  Schottky  lowering.  The  third 
term  is  responsible  for  tunneling.  A  and  b  have  been  determined  by  com¬ 
parison  with  the  experimental  results.  From  fits  to  the  experimental 

data,  a  and  b  were  found  to  be  1  x  10  e(cm^  and  2.59  x  10 

1/2 

e(cm  V)  respectively.  We  plot  in  Figure  5.4  the  barrier  heights  against 
the  oxide  fields  for  two  cases:  tunneling  lowering  effect  included  and 


Figure  5.3:  The  band  diagram  for  the  Si-SiO„  structure 
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tunneling  lowering  effect  not  included.  As  the  oxide  field  increases, 
the  barrier  height  gets  lowered  by  substantial  amounts  as  shown  in  the 
f igure . 

Surprisingly,  with  the  simple  lucky  electron  model,  Ming  was  able 
to  fit  his  experimental  results  well  for  all  the  different  devices  with 
a  consistent  set  of  parameters.  For  temperatures  at  300  K  and  77  K,  the 

O  O 

parameters  are  A  =  2.9,  X  =  91  A,  and  A  =  4.3,  X  =  108  A  respectively. 

We  discuss  in  the  following  sections  our  method  of  approach  and  look  at 
this  complicated  problem  in  more  detail. 

5 . 3  Method  of  Approach 

As  discussed  in  the  last  section,  this  is  a  transport  problem  in¬ 
volving  highly  nonuniform  spatial  electric  fields.  The  electric  field  is 
not  only  determined  by  the  static  space  charge  but  also  from  the  current 
density.  In  order  to  take  both  factors  into  account,  a  complicated  self- 
consistent  calculation,  as  discussed  in  Chapter  3,  is  required.  Fortunately, 
under  the  low  level  current  condition,  which  was  the  case  for  the  experi¬ 
ments,  the  contribution  of  the  current  is  negligible.  For  a  typical  source 

-5  -42.. 

current  of  about  5  x  10  amperes  and  a  gate  area  of  5  x  10  cm  taxen  tror. 

the  experiments,  one  finds  a  carrier  concentration  of  about  6  x  10^  cm  ^ 
if  a  saturation  drift  velocity  of  1  x  10^  cm/sec  is  assumed.  We  see  that 
the  carrier  concentration  is  orders  of  magnitude  below  the  doping  concen¬ 
tration  and  hence  its  contribution  is  totally  unimportant.  Therefore,  we 
need  only  to  include  the  spatial  variation  of  the  electric  field  due  to 
the  static  charge.  The  Monte  Carlo  approach,  including  the  nonuniform 
electric  field,  has  been  discussed  in  detail  in  Section  2.4. 
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The  total  current  coining  from  the  substrate  has  two  components: 
the  drift  component  due  to  carriers  generated  in  the  depletion  region, 
and  the  diffusion  component  due  to  the  carriers  generated  in  the  bulk 
neutral  region  and  diffusing  into  the  depletion  region.  Notice  that 
this  is  a  situation  similar  to  that  of  a  p-i-n  photodiode,  as  derived  in 
Sze's  book  [77],  The  drift  and  diffusion  component  are  expressed  as 

Jdrift  =  -^Q[  1-exp  (-axm)]  (5.5) 

and 

Jdiff  =  eVLn  exp(-axm)/(l+aLn)  (5.6) 

where  Q  is  the  inciaent  photon  flux,  x  is  the  depletion  layer  width, 

U  m 

is  the  electron  diffusion  length  and  a  is  the  absorption  coefficient. 

As  argued  by  Ning,  for  typical  values  of  between  0.02  and  0.2  cm 

-  5  -3 

(which  corresponds  to  minority-carrier  lifetimes  of  between  10  and  10 

sec),  depletion  laver  widths  x  less  than  10  cm  for  devices  used,  and 
'  m 

,  .  2  ,  3  -1 

values  ot  u  between  10  and  10  cm  for  photons  of  energy  near  the  band- 

gap,  one  finds  that  L  >  1  and  x  <<  1.  From  Equations  5.4  and  5.5,  the 

n  m 

optically  generated  current  was  completely  dominated  by  the  diffusion 
component  and  should  be  approximately  independent  of  the  depletion-layer 
width.  This  means  that  we  can  start  electrons  from  near  the  depletion 
layer  edge  to  simulate  the  diffusing-in  component  coming  from  the  bulk 
neutral  region  without  considering  the  optically  generated  carriers  inside 
the  depletion  region.  Following  the  trajectories  of  these  diffusirg-in 
electrons  with  a  Monte  Carlo  simulation,  one  shoulc  be  able  to  determine 
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Che  emission  probability  by  examining  the  percentage  of  electrons  with 
energies  greater  than  the  barrier  height  reaching  the  interface.  Un¬ 
fortunately,  this  straightforward  idea  is  too  CPU  time  intensive  be¬ 
cause  of  the  extremely  low  emission  probabilities  (10  ^  to  10  ^).  For 
an  emission  probability  of  10  ^ ,  one  has  to  simulate  a  number  much 
greater  than  10^  in  order  to  get  a  reasonable  probability  within  an  ac¬ 
ceptable  statistical  fluctuation. 

To  avoid  such  difficulties,  we  separate  our  simulation  into  two 
parts,  as  illustrated  in  Figure  5.5.  In  the  first  simulation,  as  speci¬ 
fied  by  "1"  in  the  figure,  we  start  electrons  with  a  distribution  in 
quasi  equilibrium  with  the  local  field  at  a  position  where  the  potential 
energy  is  about  twice  the  barrier  height.  This  initial  distribution  is 
automatically  generated  by  letting  the  electrons  scatter  a  few  times 
under  the  local  field  strength  before  launching  them.  Actually,  the 
electrons  get  scattered  very  often  along  the  path  and  lose  their  memories 
of  the  past  after  a  lew  scatterings.  Consequently,  if  we  start  the 
electrons  far  enough  from  the  interface,  it  does  not  matter  too  much  how 
we  start  the  initial  distribution.  Following  the  trajectories  of  the 
electrons  one  by  one,  the  energy  distribution  of  the  ensemble  is  collected 
at  position  "2"  specified  in  Figure  5.5.  We  then  start  a  second  simula¬ 
tion  for  the  electrons  in  the  high  energy  tail  of  distribution  "2"  only, 
and  collect  the  energy  distribution  of  these  electrons  at  the  interface. 
Notice  that  impact  ionizations  occur  alon,-  .,e  way,  especially  in  the 
.'egion  close  to  the  interface  where  the  electric  fielcs  are  high.  By 
calculating  the  percentage  of  electrons  in  the  high  energy  tails  as 
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illustrated  by  the  shaded  regions  in  "2"  and  "3"  of  Figure  5,5,  one  is 
then  able  to  calculate  the  emission  probability  as  follows. 

If  we  start  a  total  number  of  electrons  in  the  first  simula¬ 
tion,  with  n^  electrons  arriving  in  the  high  energy  tail  of  "2"  and 
impact  ionization  events  along  the  way,  the  probability  of  finding  an 
electron  in  the  high  energy  tail  of  "2"  can  be  written  as 


(5.  7) 


The  in  the  denominator  is  crucial  because  impact  ionization  also 
contributes  to  the  experimentally  measured  total  current.  The  same 


argument  applies  to  the  second  simulation,  and  the  probability  of  finding 
electrons  in  the  high  energy  tail  of  ”3"  can  be  written  as 


P-> 


(3.8) 


where  N,  is  the  total  number  of  electrons  used  in  the  starting  high 
energy  tail  of  "2",  m„  is  the  number  of  impact  ionization  events,  and  n., 
is  the  number  of  electrons  resulting  in  the  high  energy  tail  of  "3". 
Assuming  that  these  n,  electrons  with  energies  greater  than  the  effective 
barrier  height  (Schottkv  Lowering  and  tunnel ing  effect  included)  all  over¬ 
come  the  barrier,  the  total  emission  probability  can  be  expressed  as 

p  =  PL  *  p2  •  <5.4) 

A  large  number  of  electrons  have  to  be  considered  in  both  stages 
of  the  simulation  in  order  to  minimize  the  standard  deviation.  Depending 
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on  the  emission  probability,  an  ensemble  of  2000  to  1000  electrons  for 
the  first  simulation  and  4000  to  8000  for  the  second  simulation  is  re¬ 
quired  to  get  an  acceptable  result.  Mathematically,  the  percentage 
error  can  be  expressed  as 
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machine,  one  electron  takes  about  15  to  30  CPU  seconds  in  the  first  simu 
lation  and  5  to  10  CPU  seconds  in  the  second  simulation.  This  results 
in  about  14  to  28  CPU  hours  for  one  point! 

The  "partitioning"  of  the  real  space  is  not  unique  in  our  method. 
Care  has  to  be  exercised  in  choosing  the  boundary  conditions  for  the  two 
simulations.  In  the  first  simulation,  one  has  to  simulate  the  electrons 
all  the  way  to  the  interface  in  order  to  get  a  physically  correct  distri¬ 
bution  at  position  "2",  since  some  electrons  crossing  plane  "2"  mi  chi 
•ack  again  because  of  scatterings.  In  the  second  simulation,  iUl 
which  get  scattered  and  travel  to  the  right  siue  of  tile  so 
plane  should  be  properl,,  relaunched  according  to  tneir  : 
simulations,  electrons  hitting  the  interface  wit;.  . 
barrier  height  are  traced  a  few  more  scatter  In.-- 
sidered  to  be  collected  by  either  the  ir  ;  . 
that  as  long  as  the  boundary  cone :  •  . 

"partitioning"  does  not  uff 
the  simulation.  Tilt-  ' 

affect  the  high!  v  - 
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As  mentioned  in  the  introduction  section  of  this  chapter,  the 
purpose  for  modeling  the  emission  problem  is  twofold:  to  understand 
the  physical  details  behind  the  emission  process  and  to  determine  a 
correct  set  of  high  field  transport  parameters  to  penetrate  farther  up 
in  the  band.  Thus,  in  the  next  section,  we  concentrate  on  studying  one 
device  in  Ning's  experiment  which  gave  the  highest  emission  probability. 
This  was  the  device  numbered  as  15-2-9  and  the  potential  and  electric 
field  profiles  of  which  are  plotted  in  Figures  5.2a  and  b.  It  is  still 
too  costly  and  not  feasible  to  use  the  two  simulation  scheme  to  study 
emission  processes  of  very  low  probability. 

5.4  Results  and  Discussions 

The  three  sets  of  parameters  discussed  in  the  last  chapter  give  a 
very  different  distribution  function  in  the  high  energy  tail.  With  the 
emission  probabilities  measuring  essentially  the  characteristics  of  the 
high  energy  tail  of  the  distribution  function,  the  first  and  the  third 
sets  are  readily  eliminated.  With  too  high  a  total  scattering  rate,  the 
high  energy  tail  is  totally  suppressed,  using  the  first  set  of  parameters, 
and  rendered  the  emission  process  impossible.  With  too  low  a  total 
scattering  rate,  the  third  set  of  parameters  allows  too  high  an  emission 
probability  which  is  totally  incompatible  with  the  experiments.  In  all 
of  the  results  presented  in  this  chapter,  Set  2  parameters,  which  are  in 
between  the  two  extremes  (Keldysh's  formalism  with  P  >>  1  and  Kane’s  low 
direct  calculated  result),  were  used  in  the  simulations. 

In  Figure  5.6,  we  show  the  trajectory  of  a  typical  electron  coming 
from  the  substrate  toward  the  interface.  It  travels  along  the  band  edge 


Electron 
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Figure  5.6:  Energy  histogram  of  an  electron. 
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and  encounters  numerous  scatterings  along  the  way.  Notice  that  it  some¬ 
times  gets  back  scattered  and  travels  against  the  field.  Clearly,  this 
electron  did  not  make  it  to  the  silicon  dioxide  layer.  The  kinetic 
energy  of  the  same  electron  is  plotted  in  Figure  5.7.  The  slopes  of  the 
almost  vertical  lines  are  proportional  to  the  electric  field  at  that 
position  as  can  be  seen  from  the  equation 

AE  =  e  F  Ax  ,  (5.11) 

where  F  is  the  local  electric  field.  The  electron  sometimes  luckily 
escapes  the  scatterings  and  gains  significant  kinetic  energy.  It  loses 
energy  from  emissions  of  phonons,  traveling  against  the  field  because  of 
back  scatterings  or  suffering  impact  ionizations. 

The  calculated  emission  probability  for  device  15-2-9  is  plotted 
in  Figure  5.8  as  a  function  of  the  substrate  voltage.  The  solid  line 
gives  the  experimental  results  from  Ning  et  al.  The  open  circles  repre¬ 
sent  the  emission  probability  of  a  lowered  barrier  height  of  2.57  eV 
corresponding  to  an  oxide  field  of  2  x  10^  V/cm  as  plotted  in  Figure  5.4. 
They  are  lower  than  the  experimental  data.  In  order  to  fit  the  data, 
either  the  barrier  height  must  be  lowered  to  allow  more  electrons  to  pass 
through,  or  the  total  scattering  rate  must  be  suppressed  to  heat  up  the 
distribution  a  bit  more.  The  parameter  for  lowering  due  to  the  tunneling 
effect  has  been  obtained  from  fits  to  the  experimental  data  using  the 
lucky  electron  model.  However,  the  accuracy  of  such  a  parameter  depends 
entirely  on  how  closely  the  model  resembles  the  physical  situation.  Also, 
as  will  be  discussed  in  the  next  chapter,  the  collision  broadening  effect 
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Figure  5.7:  The  histogram  of  the  kinetic  energy  of  the 
same  electron  as  in  Figure  5.6. 
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Figure  5.8:  The  calculated  emission  probability  as  a 
function  of  the  substrate  voltage. 
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smears  out  the  sharp  edge  of  a  barrier  and  introduces  an  effective 
barrier  lowering.  If  we  allow  an  extra  .07  eV  lowering  due  to  the 
reasons  above,  we  have  the  emission  probability  represented  by  the  closed 
circles  in  Figure  5.8.  The  error  bars  attached  to  the  results  are  cal¬ 
culated  by  using  Equation  5.10.  As  discussed  at  the  end  of  the  last 
chapter,  the  computing  time  is  so  costly  that  we  cannot  afford  to  cal¬ 
culate  emission  probabilities  for  substrate  voltages  below  11  volts.  But 
if  we  extrapolate  our  calculated  data,  we  see  that  they  agree  with  the 
experimental  results. 

Figure  5.9  shows  the  energy  distributions  of  electrons  at  three 
positions  for  the  first  simulation  for  a  substrate  voltage  of  15  volts. 
The  high  energy  tail  of  the  distribution  extends  farther  up  as  we  move 
closer  to  the  interface.  In  Figure  5.10,  the  energy  distributions  for 
the  second  simulation  are  shown  for  three  positions.  Since  electrons  are 
started  with  a  high  energy  as  illustrated  in  Figure  5.5,  the  transient 
relaxation  effect  of  the  electron  energy  is  clearly  demonstrated  in  this 
figure.  As  the  electrons  move  closer  to  the  interface,  the  electron 
energy  distribution  relaxes  and  peaks  at  a  lower  value.  In  addition,  the 
high  energy  tails  broaden  due  to  the  acceleration  of  high  fields  close  to 
the  interface.  About  half  of  the  electrons  in  the  high  energy  tail  are 
in  the  second  conduction  band. 

To  look  at  the  effect  of  impact  ionization,  we  plot  in  Figure  5.11 
the  average  ionization  coefficient  as  a  function  of  the  surface  electric 
fields.  The  dashed  line  represents  the  steady  state  experimental  re¬ 
sults  from  Lee  et  al.  [58]  and  the  triangles  are  the  calculated  results 
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Electron  Energy  (eV) 


:  The  energy  distribution  of  electrons  at  three  positions 
for  the.  first  simulation. 
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Figure  5.10:  The  energy  distribution  of  electrons  at  three 
positions  for  the  second  simulation. 
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Figure  5.11:  The  electron  ionization  coefficient  plotted 


against  the  surface  field. 
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corresponding,  from  left  to  right,  to  substrate  voltages  of  19,  17,  15, 
and  13  volts.  The  electric  field  is  highest  at  the  interface  and  de¬ 
creases,  rapidly  moving  toward  the  substrate.  The  average  electric 
field  the  electrons  undergo  is  much  smaller  than  the  maximum  surface 
field.  With  the  triangles  lying  only  slightly  lower  than  the  steady  state 
curve,  it  seems  to  suggest  that  the  transient  ionization  process  under  an 
increasing  field  condition  is  enhanced.  In  Figure  5.12,  the  scattering 
mean  free  path  in  the  field  direction  is  plotted  against  the  substrate 
voltages.  We  see  that  the  scattering  mean  free  path  is  not  a  constant 
but  a  strong  function  of  the  substrate  voltage  (or  the  electric  field). 

O 

They  are  much  shorter  than  the  constant  mean  free  path  of  91  A  used  by 
Ning  et  al. 

To  understand  this  discrepancy,  let  us  reexamine  the  idea  of  the 

lucky  electron  theory.  If  momentum  space  is  considered,  then  the  idea 

of  the  lucky  electron  theory  is  to  find  the  probability  of  an  electron 

— ► 

starting  from  a  state  k^,  not  encountering  any  scattering  in  a  time  in¬ 
terval  t,  and  ending  up  in  a  state  k  given  as 
-► 

k  =  •  (5-12) 

The  idea  is  conceptually  simple  and  mathematically  easily  manageable. 

But  the  most  important  and  difficult  part,  which  is  always  neglected, 
is  the  initial  distribution  of  kQ.  In  the  electron  emission  model,  what 
has  been  neglected  is  the  initial  distribution  at  the  starting  position 
where  the  potential  energy  equals  the  emission  barrier  as  shown  in 
Figure  5.3.  There  the  electrons  should  have  an  energy  distribution  like 
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The  scattering  mean  free  path  in  the  field  direct 
plotted  as  a  function  of  substrate  voltage  for 
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those  shown  in  Figure  5.9  with  an  average  energy  of  about  1  eV.  If  an 
average  phonon  energy  of  0.05  eV  is  considered,  then  the  electrons  can, 
on  the  average,  suffer  20  phonon  emissions  and  still  make  it  to  the 
silicon  dioxide  layer.  Although  King's  lucky  electron  model  agrees 
with  the  experimental  data,  we  doubt  that  one  can  relate  direct  phys¬ 
ical  meanings  to  the  essentially  adjustable  parameters.  Mathematically, 
however,  his  model  is  simple  and  appealing  and  can  be  directly  applied 
to  the  modeling  of  devices. 

5. 5  Conclusion 

We  have  shown  that  a  transient  Monte  Carlo  simulation  can  be  used 
to  correctly  describe  the  emission  process  of  electrons  into  the  silicon 
dioxide  layer.  This  Monte  Carlo  method  takes  into  account  the  non- 
uniform  spatial  electric  field  which  is  the  most  important  characteristic 
of  the  system.  In  this  study,  the  inclusion  of  a  realistic  band  structure 
is  essential  because  this  transport  process  involves  very  high  energy 
electrons  not  describable  by  an  effective  mass  approximation.  Also,  the 
second  conduction  band  plays  an  important  role  since  about  half  of  the 
emissions  are  from  this  band.  We  are  able  to  determine  the  right  set  of 
high  field  transport  parameters  from  this  study.  This  should  be  con¬ 
sidered  as  a  major  success  which  improves  the  understanding  of  high  field 
transport  processes,  like  impact  ionization  phenomenon,  and  enables  us  to 
penetrate  farther  up  into  bands.  Previous  theories  for  impact  ionization 
using  Keldysh's  formalism  with  P  >>  1  have  shown  to  be  incorrect  in  that 
the  high  energy  tail  of  the  distribution  function  is  overly  suppressed. 

The  subtle  balance  between  the  high  energy  phonon  scattering  rate  and 
impact  ionization  scattering  rate  as  described  in  Chapter  4  is  by  no  means 
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optimized  because  of  the  costly  computing  time  which  limits  our  freedom 
to  "play  around"  with  the  parameters.  Unlike  previous  models  for  the 
electron  emission,  impact  ionization  is  naturally  included  in  our  cal¬ 
culation  . 

We  have  also  shown  that  the  physics  behind  the  emission  process 
can  not  be  explained  by  the  simple  lucky  electron  model.  Although  the 
lucky  electron  model  consistently  fits  the  experimental  results,  one 
could,  at  best,  look  at  it  as  a  curve  fitting  scheme.  For  device  modeling 
purposes,  this  simple  model  is  certainly  the  best  one  to  use  mathematically. 
In  Chapter  6,  we  discuss  a  better  way  to  calculate  the  high  energy  phonon 
scattering  rate,  including  the  quasi-particle  self-energy  effect. 
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CHAPTER  6 

QUANTUM  TRANSPORT  IN  SEMICONDUCTORS 


6. 1  Introduction 

The  theoretical  basis  for  the  transport  study  has  been  the  semi- 
classical  Boltzmann  equation.  This  equation  is  conceptually  simple  and 
mathematically  manageable.  On  a  large  time  and  length  scale,  classical 
descriptions  of  the  motion  of  the  particles  are  usually  acceptable. 

But  as  the  length  scale  of  devices  shrinks  to  the  submicron  range,  where 
extremely  high  electric  fields  are  unavoidable,  one  has  to  question  the 
validity  of  the  semiclassical  Boltzmann  equation. 

The  formalism  of  the  quantum  transport  theory  was  basically  de¬ 
veloped  in  the  pioneering  work  of  Kohn  and  Luttinger  (1957,  1958) [  78,79]  , 
Kubo  (1957) [80],  Dresden  (1961) [81] ,  Chester  (1963) [82] ,  Kubo  (1966) [83], 
and  Luttinger  (1968) [84].  The  formalism  is  quantum  mechanically  rigorous, 
but  unfortunately,  it  is  neither  conceptually  simple  nor  mathematically 
manageable.  The  interesting  work  of  Scott  and  Moore  in  1972  [85]  showed 
that  the  quantum  transport  equation  can  be  reduced  to  a  quasi-particle 
Boltzmann  transport  equation.  This  is  encouraging  in  that  one  needs  only 
to  include  in  the  Boltzmann  transport  theory  the  quasi-particle  character 
which  amounts  essentially  to  a  self-energy  renormalization  effect.  In 
1973,  Barker  rederived  the  quantum  transport  theory  for  high  field 
transport  using  a  superoperator  technique  for  semiconductors,  and  dis¬ 
cussed  in  his  work  the  possible  first  order  quantum  corrections  to  amend 
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the  breakdown  of  the  Boltzmann  equation  [86].  The  results  of  his  work 
also  provide  criteria  for  using  the  Boltzmann  equation. 

This  chapter  is  devoted  to  the  discussion  of  the  classical  and  quantum 
transport  equations  and  setting  up  confidence  limits  for  the  application 
of  the  Boltzmann  equation.  In  Section  6.2,  the  assumptions  of  the 
Boltzmann  equations  are  discussed  and  the  usual  approximations  made  in 
solving  the  equation  are  examined.  In  Section  6.3,  the  basis  of  the  quan¬ 
tum  transport  equation  is  introduced  and  various  quantum  effects  due  to 
the  high  electric  field  and  the  strong  coupling  of  the  electron-phonon 
system  are  discussed.  The  quasi-particle  concept  is  introduced  and  its 
self-energy  due  to  the  interaction  with  phonons  is  calculated.  We  also 
show  how  the  scattering  rate  is  obtained  from  the  self-energy.  A  quantum 
Monte  Carlo  scheme  is  proposed.  In  Section  6.4,  we  show  the  transport 
results  with  the  quantum  corrections  and  assess  the  importance  of  the 
various  quantum  effects. 


6. 2  Assumptions  of  the  Boltzmann  Equation 

The  semiclassical  Boltzmann  equation  is  of  the  form 
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(1)  Electrons  are  independent  classical  particles  so  that  a  distri- 
bution  function  f(r,k;t)  can  be  defined  in  classical  phase  space 
(r,k) . 

(2)  The  electronic  states  are  stationary  (long  life  time)  with  a  well- 
defined  momentum  k. 

(3)  The  effects  of  impurities  and  lattice  vibrations  (phonons)  can  be 
considered  as  perturbations  causing  weak  scattering  among  the 


Bloch  states. 
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(A)  The  external  field  solely  accelerates  the  electrons  between  col¬ 
lisions  and  has  no  effects  on  the  electronic  states  and  the 
scattering  events. 

(5)  The  scattering  events  are  considered  to  be  local  both  in  time  and 
space.  The  duration  of  a  collision  is  negligible.  Thus,  the  tran¬ 
sport  processes  are  viewed  on  a  coarse-grained  time  scale  x  >>  x  , 

where  x  is  the  collision  duration, 
c 

From  the  assumptions,  we  see  that  the  basic  objection  in  the  use 
of  the  Boltzmann  equation  is  the  exclusion  of  the  uncertainty  principle 
in  the  formalism.  This  raises  the  question  whether  the  distribution 
function  f(r,k;t)  can  always  be  meaningfully  defined  in  the  classical 
phase  space.  If  a  steady  state  solution  is  of  interest,  then  f(k)  can 
be  properly  defined  without  violating  the  position-momentum  uncertainty 
principle.  But  there  is  always  the  energy-time  uncertainty.  Care  has 
to  be  exercised  in  using  the  latter  uncertainty  relation,  since  time  is 
not  an  operator  [87].  Let  us  follow  the  qualitative  discussion  of 
Dresden  [81]  to  explain  how  the  uncertainty  principle  affects  the  concept 
of  the  distribution  function. 

For  a  given  distribution  function  f(E)  and  a  state  density  function 
p(E),  the  number  density  of  electrons  is  given  as 

n(E)  dE  =  f (E)  p(E)  dE  .  (6.7) 

n(E)  is  interpreted  as  the  number  of  electrons  in  an  interval  dE  around 
energy  E.  The  number  density  changes  rapidly  in  a  region  kT  around  the 
average  energy  of  electrons  in  a  semiconductor.  (In  metals,  it  is  the 
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Fermi  energy. )  In  order  to  meaningfully  define  the  distribution  func¬ 
tion,  one  should  be  able  to  define  the  number  of  electrons  in  energy 
ranges  small  compared  to  kT.  The  basic  uncertainty  in  the  energy  of  the 
electrons  is  of  the  order  of  -fi/x,  where  1/x  is  the  total  transition  rate. 
Thus,  one  must  demand  that 

-  «  kT  .  (6.8) 

Peierls  has  derived  this  relation  quantitatively  from  time  dependent 
perturbation  theory  [88]  but  he  has  also  shown  that  this  relation  can  be 
relaxed.  In  Equation  6.8,  the  relevant  energy  should  be  E,  the  electron 
energy,  instead  of  the  thermal  energy,  kT.  Thus,  it  is  only  required 
that 


Y  «  E  .  (6.9) 

Landau  has  also  shown  that  for  elastic  scatterings  only,  the  Boltzmann 
equation  can  be  derived  without  the  assumption  of  Equation  6.7.  Qualita¬ 
tively,  Equation  6.8  reads  that  the  uncertainty  in  energy  should  be 
negligible  compared  to  the  total  kinetic  energy  of  the  electrons.  This 
actually  holds  true  for  most  of  the  semiconductors. 

Much  of  the  incorrect  nature  of  the  classical  Boltzmann  equation  can 
be  removed  by  introducing  the  notion  of  a  quasi  particle,  which  will  be 
discussed  in  the  following  sections.  The  more  serious  problems  actually 
lie  in  the  usual  approximations  for  solving  the  Boltzmann  equation.  For 
transport  problems  in  semiconductors,  these  are  the  effective  mass  approx¬ 
imation  (good  for  low  field  but  not  acceptable  for  high  field  transport). 
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the  neglect  of  electron-electron  interactions  (good  for  low  electron 
density  only) ,  the  use  of  the  Golden  Rule  and  the  Born  approximation  to 
calculate  the  scattering  rate  (good  for  low  scattering  rate  only),  and 
the  neglect  of  the  effect  of  electric  fields  on  scatterings. 

As  will  be  discussed  in  the  following  sections,  it  turns  out  that 
the  self-energy  effect  of  a  quasi  particle  plays  the  most  important  role 
in  modifying  the  classical  Boltzmann  transport  equation. 

6.3  Quantum  Transport  and  Quasi  Particles 

6.3.1  Introduction 

The  basis  of  the  quantum  transport  theory  is  the  quantum  Liouville 
equation : 

iff  •  ^  =  [H,n]  ,  (6.10) 

where  0  is  the  density  matrix  of  the  system  including  both  the  carrier 
part  and  the  phonon  part,  and  H  is  the  full  Hamiltonian  including  the 
external  force.  For  small  external  forces,  the  linear  response  theory 
developed  by  Kubo  [83]  can  be  applied  to  solve  Equation  6.10.  In  the  non¬ 
linear  response  regime  caused  by  large  external  forces.  Equation  6.10 
becomes  a  mathematical  nightmare  in  that  the  linearization  of  the  equation 
can  not  be  made  and  the  phonon  part  and  the  carrier  part  of  the  density 
matrix  can  not  be  easily  separated.  The  superoperator  technique  developed 
by  Zwanzig  [89]  serves  as  a  neat  mathematical  tool  for  finding  the  formal 
solution  of  Equation  6.10.  But  in  order  to  get  from  the  formal  solution 
a  mathematically  manageable  transport  equation,  some  assumptions, 
not  immediately  justifiable,  are  made  in  the  derivation.  The  only  true 
justification  of  the  assumptions  made  lies  in  the  fact  that  the  Boltzmann 
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equation  can  be  reconstructed  from  the  quantum  transport  equation  in  the 
classical  limit. 

Aside  from  the  purely  theoretical  interest  in  developing  the  quantum 
transport  theory,  the  purpose  has  been  to  look  at  possible  first  order 
correction  terms  to  amend  the  classical  Boltzmann  equation.  Scott  and 
Moore  [85]  did  an  interesting  study  on  the  quantum  electron- impurity 
transport  problem.  They  used  a  Green's  function  technique  and  formally 
derived  a  quasi-particle  Boltzmann  equation.  Their  result  suggests  that 
the  most  important  quantum  correction  in  the  transport  problem  is  to 
bring  in  the  notion  of  a  quasi-particle  for  the  Boltzmann  equation. 

6.3.2  Summary  of  Barker's  Results  [86] 

Barker  has  adopted  the  superoperator  technique  and  rederived  the 
quantum  transport  theory  for  high  field  transport  in  semiconductors.  He 
has  investigated  the  effect  of  a  high  electric  field  on  the  scattering  of 
the  carrier.  Here,  we  summarize  his  results  and  discuss  qualitatively 
their  physical  implications. 

( i)  Self-Energy  Effect 

As  discussed  in  Section  6.2,  the  scattering  processes  are  treated  as 
real  transitions  between  sharp,  unperturbed  momentum  states  in  the  classical 
Boltzmann  transport  picture.  The  very  existence  of  scattering  and  the 
presence  of  many  scatterers  make  the  assumption  of  an  isolated  scattering 
unrealistic.  In  fact,  the  carrier  propagates  in  a  perturbed  state  which 
is  controlled  by  virtual  scattering  processes.  The  net  effect  is  that  the 
energy  of  a  state  is  shifted  by  an  amount  which  is  the  real  part  of  the 
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self-energy,  and  Che  life  time  of  a  state  is  finite  and  is  proportional 
to  the  imaginary  part  of  the  self-energy.  The  state  is  collision 
broadened.  Because  of  the  self-energy  effect,  the  energy  conserving 
delta  function  becomes  a  smeared-out  Lorentzian 

5(E(k)  -  E(k-q)  -  «*+)  -  ±  — - — - r(K’^T-2  ~ * 

q  (E (k)  -  E(k-q)  -  hur>-A(k,q))Z  +  T  (k,q) 

q 

(6.11) 


where 

A(k,q)  =  A(k,E(k-q) )  -  A(k-q,E(k) ) ,  (6.12) 

and 

"(k,q)  *  F  (k,E (k-q) )  +  F(k-q,E(k)).  (6.13) 

The  quantity  A  is  the  real  part  of  the  self-energy  and  F  is  the  imaginary 
part  of  the  self-energy.  Notice  that  when  A  and  F  approach  zero,  which 
is  the  case  for  weak  coupling,  the  Lorentzian  structure  goes  back  to  the 
form  of  a  delta  function.  The  physical  implication  of  the  self-energy  ef¬ 
fect  is  that  the  particle  is  scattered  to  a  range  of  quasi  particle  states 
according  to  the  Lorentzian  distribution. 

What  has  been  described  is  essentially  the  characteristics  of  a  quasi¬ 
particle.  We  will  come  back  to  more  detail  when  we  introduce  the  quasi¬ 
particle  concept. 

(ii)  Intra-Coll lsional  Field  Effect 

In  an  actual  scattering  event,  the  interaction  between  the  particles 
takes  a  finite  time.  In  the  presence  of  an  electric  field,  the  carrier 
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gains  (or  looses)  energy  during  the  finite  collision  duration.  If  the 
energy  gained  (or  lost)  during  the  collision  duration  is  much  greater 
than  the  collision  broadening  of  the  state  of  the  carrier,  then  modifi¬ 
cations  need  be  made  to  the  scattering  rate.  The  energy  gained  (or  lost) 
during  the  collision  is  expressed  as 

E  =  | eF  •  v I  •  X  ,  (6.14) 

cd  1  1  c 


where  F  is  the  electric  field,  v  is  the  velocity  of  the  carrier,  and 

is  the  collision  duration.  For  a  rough  estimate,  the  collision  dura¬ 
tion  can  be  written  as 


(6.15) 


where  v  is  again  the  velocity  of  the  carrier,  and  A  , 
length  of  the  carrier,  is  given  as 


the  deBroglie  wave 


(6.16) 


where  k  is  the  crystal  momentum  of  the  carrier.  In  other  words,  the  collision 
duration  is  just  one  period  of  the  particle  wave.  Quantum  mechanically, 
one  can  think  of  the  collision  duration  as  the  time  it  takes  for  a  par¬ 
ticle  to  fully  manifest  its  wave  characteristics.  Mathematically,  the 
collision  broadened  scattering  rate,  calculated  by  using  Equation  6.11,  is 
essentially  unmodified  by  the  presence  of  an  electric  field  if  the  follow¬ 
ing  criterion  is  met: 


r 


<  i  , 


(6.17) 


Ill 


where  E  .  is  given  by  Equation  6.14  and  f  is  given  bv  Equation  6.13. 
cd 

In  the  weak  coupling  regime  (low  carrier  energy,  low  scattering 
rate)  where  T  is  small,  the  criterion  given  by  Equation  6.17  might  be 
violated.  But  the  carrier  distribution  is  heated  up  significantly  under 
high  electric  fields  so  that  most  of  the  carriers  reside  at  high  energies, 
and  hence  in  the  strong  scattering  regime,  where  T  is  very  large.  It  is 
unlikely  that  Equation  6.17  is  violated  for  carriers  in  the  high  energy 
range.  In  Section  6.4,  we  present  the  Monte  Carlo  result  to  assess 
quantitatively  the  intra-collisional  field  effect. 


6.3.3  Quasi  Particle 

For  a  noninteracting  system,  the  wave  function  of  a  particle  in  a 
momentum  state  k  evolves  as 


(6.18) 


where  E(k)  is  the  energy  of  the  free-particle  state  k.  The  propagation 
of  this  particle  is  undamped  and  the  frequency  of  oscillation  is  just 
the  free-particle  energy  E(k).  If  we  take  the  Fourier  transform  of 
Equation  6.18,  we  obtain  the  spectral  density  function  as 

— ► 

A  (to ,  E )  ~  6(( o  -  ^4r^-)  .  (6.19) 

For  an  interacting  system,  the  energy  of  the  particle  in  a  momentum 
state  k  is  modified.  The  interaction  introduces  an  extra  energy  to  the 
particle  and  hence  the  evolution  of  the  wave  function  of  the  particle  is 


determined  by 
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_ .  (E (k)  +  I(k))  t 

4»(t)  ~  e  *  (6.20) 

where  I(k)  is  called  Che  self-energy  of  Che  parCicle  and  is  in  complex 
form  expressed  as 

ZT(k)  =  A(k)  -  ixr (k)  .  (6.21) 

Thus,  Che  wave  function  evolves  as 

1  E(k)  -6  A  (k)  _T(k) 

Ut)~e~  *  e  *  .  (6.22) 

The  real  part  of  the  self-energy  introduces  a  level  shift  of  the  eigen¬ 
states,  and  the  imaginary  part,  which  is  related  to  the  level  broadening, 
damps  the  magnitude  of  the  wave  function.  The  magnitude  square  of  the 
wave  function,  which  represents  the  probability  of  a  particle  in  a  state, 
follows 

_ t_ 

|iMt)|2  ~  e  \  ,  (6.23) 

where 

2F  (k) 

is  the  life  time  of  a  particle  in  the  momentum  state  k.  In  first 

order  approximation,  1/t^  corresponds  to  the  scattering  rate  of  the 

— ► 

particle  in  the  state  k.  This  will  be  justified  when  we  formally  derive 
the  self-energy.  If  we  take  the  Fourier  tra  .lot  of  Equation  6.22  and 
change  sign  in  the  second  exponent  for  t  <  0,  we  have 
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I 
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One  sees  that  the  spectral  density  function  takes  the  Lorentzian  fom. 
and  peaks  at  the  energy  E(k)  +  A.  The  state  is  broadened. 

We  see  that  in  an  interacting  system,  the  state  has  a  natural  life 
time  and  the  particle  decays  out  of  its  initial  state.  In  order  to 
still  use  the  concept  of  a  particle,  one  defines  the  concept  of  a  quasi 
particle.  A  quasi  particle  in  a  momentum  state  k  is  defined  as  a  par- 

— ►  —V 

tide  with  momentum  k,  energy  E(k)  +  A  and  a  broadened  spectral  density 
function  given  by  Equation  6.25.  The  concept  is  justified  for  times  t 
which  satisfy 


« 

E 


< 


(6.26) 


What  has  been  described  are  simple,  qualitative  concepts  of  a  quasi 
particle.  For  a  rigorous  discussion,  the  reader  is  referred  to  Pine’s 
book,  The  Many  Body  Problem  [90].  The  important  task  remaining  is  to 
find  the  self-energy  which  characterizes  the  properties  of  a  quasi 
particle. 

We  start  with  an  approximation  to  the  self-energy  E  which  neglects 
the  vertex  correction  [91].  It  is  represented  by  the  diagram  shown  in 
Figure  6.1.  The  double  line  represents  the  full  propagator  which  is  the 
dressed  Green's  function  of  the  electron.  It  emits  phonons  and  reabsorbs 


them  back  by  virtual  processes.  The  self-energy  of  the  electron  due  to 
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which  can  be  evaluated  by  contour  integration.  The  integrand  has  three 

poles  in  the  complex  plane.  The  two  poles  fTii)  =  -  -ftorS+id  and 

•R u;  =  E-E(k-q)  -  Rel  +  i(6-lml)  are  on  the  upper  half  plane  (F  =  -ImZ, 

which  is  proportional  to  the  life  time,  is  always  positive),  and  the 

pole  fiio  =  -ttoJ  -  i5  is  on  the  lower  half  plane.  One  can  either  close  the 

q 

contour  in  the  upper  or  lower  half  planes,  depending  on  which  is  more  con¬ 
venient.  By  closing  the  contour  in  the  lower  half  plane,  one  finds 

I  =  - —  -  (6.31) 

E  -  -ftur'-  -  E(k-q)  -  I  (k-q  ,  E-Rur*)  +  i6 

q  q 

bv  applying  the  residue  theorem  [95] .  The  self-energy  expression  is 
thus  simplified  to 


:(k,E)  = 


a  (k+q') 


•I  ( 2tt ) 3  E-ftu^-E(q')  -  Z (q '  ,E-fiw°  )  +  io 
q  k+q’ 


(6.32) 


where  q'  =  k-q  and  g  (q)  lumps  all  the  factors  together.  Let  us  consider 

the  case  of  the  optical  phonon  scattering  with  and  assume  that 

g (q )  is  a  constant  instead  of  a  function  of  q,  i.e.,  the  usual  approxi- 

— ► 

mation  of  a  constant  deformation  potential.  One  sees  that  Z(k,E)  be- 
comes  a  function  of  the  energy  only  after  the  integration  of  q' . 

Utilizing  the  relation 


r  3  _  r 

-  -  q  F (E (q) )  =  I  dC  p  (E)  F (E)  , 

1  ( 2tt)  J  1 


(6.33) 
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where  F  is  any  function,  p(E)  is  the  density  of  states  and  the  crystal 
volume  is  taken  as  one,  Equation  6.32  then  takes  the  form 


£ (E)  =  g2 


_ dE 1  p  (E* ) _ 

E-tfw  -  E'  -  Z(E  -  )  +  id 

op  op 


(6.34) 


Equation  6.34  can  be  solved  numerically  using  the  density  of  states 
calculated  from  the  pseudopotential  band  structure  [37], 

Let  us  now  consider  the  weak  coupling  limit  (g  <<  1)  where  £  in 
the  denominator  can  be  neglected,  i.e.,  an  equation  for  the  lowest 
order  self-energy.  By  use  of  Dirac’s  formula  [95] 

- - -  =  P  ~~  ±  iTrS(x'-x)  ,  (6.35) 

x'-x  ±  id  x  “x 

Equation  6.34  becomes 


1(E) 


’  dE'g2?(E’) 

E-tkj  -E' 

Op 


iTr g  p(E-hw  ) 
op 


(6. 36) 


For  Z(E)  =  A(E)  -  ir(E), 


A(E)  =  P 


dE,g^p(E' ) 

E— Hoj  -E' 
op 


(6.37) 


and 


r(E)  =  ng2  p(E-ftw  )  . 

op 


(6.38) 
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Note  that  the  real  part  A(E)  is  just  the  energy  correction  of  the  second 
order  perturbation  theory  [96]  and  the  imaginary  part  T(E)  is  similar  to 
the  form  of  the  Golden  rule  in  the  Born  approximation.  As  discussed 
earlier  in  the  text,  2h/T(E)  is  the  life  time  of  the  state.  The  self- 
energy  I(E(k)),  which  is  due  to  all  possible  virtual  scatterings  of  the 
same  initial  and  final  state,  represents  the  forward  scattering  amplitude. 
In  the  approximation  that  the  vertex  part  of  the  diagram  can  be  neglected 
[97],  Z(E(k))  is  simply  <k  |  T  J  k>,  where  T  is  the  transition  matrix  [98]. 

If  we  assume  the  scattering  is  elastic,  then  the  scattering  probability 
is  conserved  and  the  optical  theorem  [98] 


tot  hv, 


Im  <kj  T  k> 


hv, 


Im  Z (E (k) ) 


can  be  applied  to  give  the  total  scattering  rate  as 


(6.39) 


1 


tot 


2F (E (k) ) 
h 


(6.40) 


One  sees  that  the  scattering  rate  is  just  the  inverse  life  time  of  the 
state.  Thus,  the  total  scattering  rate  is  given  as 


7717  =  'hL  g2  p(E-hwop)  ’  (6-41) 

D 

which  is  exactly  the  Golden  rule  in  the  Born  approximation  for  phonon 
2  2 

emission  with  Ih'|  =  g  ,  the  squared  matrix  element  as  given  in  Appendix  2. 
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Since  we  start  out  with  the  ground  state  Green’s  function,  only  phonon 
emission  can  occur  in  the  scattering  process. 

What  has  been  derived  is  valid  for  high  energy  scattering  processes 
where  the  zero-point-lattice  serves  as  a  good  approximation.  By  solving 
Equation  6.34,  the  total  scattering  cross  section  can  be  calculated  from 
Equation  6.39  in  the  elastic  limit.  For  impurity  scatterings.  Equation 
6.39  is  exact.  In  a  real  solid,  there  is  more  than  one  scattering 
mechanism,  for  example,  the  scatterings  of  different  phonon  types,  the 
impact  ionization  scattering,  etc.  The  self-energy  formalism  discussed 
above  is  only  true  for  one  type  of  scattering.  One  can  treat  separately 
the  different  scattering  mechanisms  and  sum  algebraically  the  self¬ 
energies,  assuming  that  the  cross  diagrams  for  different  scattering  types 
can  be  neglected.  The  justification  for  the  assumption  is  an  open  ques¬ 
tion  which  needs  to  be  further  pursued.  We  propose  in  the  next  section 
a  quantum  Monte  Carlo  procedure  in  the  spirit  of  a  quasi-particle  picture 

6.3.4  Quantum  Monte  Carlo 

Up  to  now,  we  have  discussed  that  the  most  important  quantum  cor¬ 
rection  to  the  transport  problem  is  to  bring  in  the  notion  of  a  quasi¬ 
particle.  As  mentioned  in  Chapter  2,  it  has  been  proven  that  a  Monte 
Carlo  simulation  solves  the  Boltzmann  equation  [16,20].  With  the  quasi¬ 
particle  Boltzmann  equation  being  formally  the  same  as  the  classical 
Boltzmann  equation,  it  is  also  solved  by  the  Monte  Carlo  simulation  pro¬ 
vided  the  quasi  particle  characteristics  are  properly  included  in  the 
simulation.  The  self-energy,  which  is  the  main  characteristic  of  a  quasi 
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particle,  has  two  effects  on  a  momentum  state:  it  shifts  the  energy  of 
the  state  by  an  amount  equal  to  its  real  part,  and  it  collision-broadens 
the  state  with  a  Lorentzian  structure  whose  half  width  equals  its  im¬ 
aginary  part.  For  our  Monte  Carlo  method,  including  a  realistic  band 
structure,  the  contribution  from  the  real  part  can  be  neglected,  since  it 
is  partly  accounted  for  in  the  empirical  pseudopotential  band  structure. 

As  for  the  imaginary  part,  the  broadened  spectral  density  function 
(Equation  6.25)  results  in  a  smeared-out  Lorentzian  structure  (Equation 
6.11).  This  takes  the  place  of  the  energy  conserving  delta  function  in 
the  scattering  process  as  derived  by  Barker  [86] .  This  factor  can  be 
taken  care  of  by  scattering  the  electron  into  a  range  of  quasi-particle 
final  states  according  to  the  Lorentzian  structure.  Assuming  a  constant 
coupling  as  discussed  in  the  last  section,  the  self-energy  is  a  function 
of  the  energy  only.  Thus,  the  joint  real  and  imaginary  parts  are  rewritten 
as 

A  =  A (E (k-q) )  -  A(E(k) )  (6.42) 

and 

f  =  r (E (k-q) )  +  F(E(k))  .  (6.43) 

We  see  that  the  joint  broadening  depends  both  on  the  energy  of  the  initial 
and  the  final  states.  It  is  very  difficult  to  realize  this  distribution 
in  the  Monte  Carlo  simulation  effectively  and  efficiently.  Our  simplify¬ 
ing  procedure  is  that  we  consider  the  initial  state  broadening  only,  and 
scatter  the  electron  into  an  energy  range  of  2I\,  the  full  width,  according 


121 


to  the  distribution.  The  partial  justification  for  scattering  only  into 
an  energy  range  of  the  full  width  is  because  the  Lorentzian  tail  extends 
both  far  up  and  down  in  energy  where  scattering  mechanisms  other  than 
the  one  responsible  for  this  broadening  might  be  important,  and  this 
formula,  derived  for  a  single  scattering  mechanism,  might  not  be  valid. 

The  equation  of  motion  for  the  particle  stays  the  same,  as  has  been  shown 
in  [85].  So  the  free  flight  time  can  be  found  the  same  way  by  solving 
Equation  2.1,  except  that  one  uses  the  scattering  rate  calculated  from 
the  self-energy  in  the  high  energy,  strong  coupling  regime. 

The  above  proposed  scheme  is  by  no  means  complete.  It  attempts  to 
serve  as  a  first  step  toward  the  further  understanding  of  the  quantum 
transport  problem. 

6.4  Results  and  Discussion 

We  discuss  first  the  results  of  numerical  solutions  to  Equation  6.34 
for  the  self-energy,  using  an  empirical  pseudopotential  band  structure 
[37,99]. 

In  Figure  6.2,  we  show  the  result  of  the  self-energy  due  to  the 
intervalley  electron-phonon  interaction  in  GaAs.  In  this  calculation,  the 
intervalley  phonon  energy  is  assumed  to  be  28  meV,  an  average  of  the 

O 

different  intervalley  phonon  energies,  and  the  coupling  constant  g  =  0.06 
is  selected  such  that  the  low  energy  scattering  rate  agrees  with  what 

describes  the  Gunn  effect  [25,36].  The  solid  line  represents  the  shifted 

2 

density  of  states  pCE-ftw^)  multiplied  by  ng  ,  which  is  the  imaginary  part 
of  the  lowest  order  self-energy  as  described  by  Equation  6.38.  The  dashed 
line  represents  the  imaginary  part  of  the  full  order  self-energy  as 


Scattering  Rate  (eV) 


Energy  E(eV) 


Figure  6.2:  The  scattering  rates  (the  imaginary  part  of  the  self¬ 
energy)  plotted  as  a  function  of  the  electron  energy. 

The  solid  curve  corresponds  to  the  first  order  self¬ 
energy  and  the  dashed  curve  corresponds  to  the  full  order 
•*lf”*nergy.  The  dash-dotted  curves  are  that  used  in  the 
calculations.  The  inset  shows  the  real  part  of  the  self¬ 
energy. 
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described  by  Equation  6.34.  We  see  that  the  full  order  self-energy  cor¬ 
rection  is  negligible  for  energies  below  ~0.6  eV.  It  is  broadened  and 
suppressed  by  about  20%  around  the  peak  of  the  density  of  states.  The 
real  part  of  the  full  order  self-energy,  which  is  partly  taken  care  of 
in  the  empirical  pseudopotential  calculation,  is  shown  in  the  inset  of 
the  figure.  Also  included  in  the  figure  are  the  scattering  rates  used 
in  the  previous  Monte  Carlo  simulation  by  Shichijo  et  al.  (high  dash- 
dotted  curve)  and  the  present  Monte  Carlo  simulation  (low  dash-dotted 
curve)  which  will  be  discussed  later.  To  illustrate  the  coupling- 
constant  dependence  of  the  scattering  rate,  we  plot  in  Figure  6.3  the 
imaginary  part  of  the  self-energy  F  for  five  coupling  constants.  We 

find  that  the  scattering  rate  at  the  density-of-states  peak  grows  al- 

2 

most  linearly  with  g  within  the  selected  range. 

There  are  many  different  phonon  types  for  the  intervallev  scattering 

in  silicon  [18].  As  discussed  in  Section  6.3.3,  it  is  not  clear  how  one 

includes  the  contribution  of  the  mixed  interaction  with  different  phonon 

types  to  the  self-energy.  Therefore,  we  pick  a  phonon  energy  of  65  meV 

2 

and  a  coupling  constant  of  g  =  .07  to  look  qualitatively  at  the  self¬ 
energy  effect.  Figure  6.4  shows  the  imaginary  part  of  the  self-energy 
for  two  conduction  bands  in  Si.  The  solid  and  the  dashed  curves  in  this 
figure  correspond  to  that  described  for  Figure  6.2.  The  scattering  rate 
from  the  full  self-energy  is  again  broadened  and  suppressed  around  the 
two  peaks.  Note  that  it  is  suppressed  by  over  50%  at  the  second  peak. 

The  real  part  of  the  self-energy  is  shown  in  the  inset  of  the  figure. 
Qualitatively,  our  empirical  scattering  rate  shown  in  Figure  4.9  agrees 
with  the  self-energy  scattering  rate  for  energies  below  3  eV. 
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Figure  6.4:  The  scattering  rate  (the  imaginary  part  of  the  self-energy) 
plotted  as  a  function  of  energy.  The  solid  line  is  the 
first  order  self-energy  (the  Golden  rule  in  the  Born  approx¬ 
imation)  and  the  dashed  line  is  the  full  order  self-energy. 
The  inset  shows  the  real  part  of  the  self-energy. 
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The  coupling-constant  dependence  of  the  scattering  is  shown  in  Figure  6.5. 
The  qualitative  behavior  is  the  same  as  that  of  Figure  6.3  for  GaAs . 

In  our  quantum  Monte  Carlo  simulation,  we  still  adopt  the  empir¬ 
ical  scattering  rate  as  discussed  in  Chapters  4  and  5  for  Si  and  the 
lower  dash-dotted  curve  shown  in  Figure  6.2  which  fits  the  experimental 
results  for  GaAs  better.  In  Figure  6.6,  the  impact  ionization  coefficient 
a  for  GaAs  is  shown  for  both  the  quantum  and  the  classical  Monte  Carlo 
simulations.  All  results  fall  into  the  cross-hatched  region  which  in¬ 
dicates  the  range  of  the  available  experimental  data  [100-103].  We  see 
that  the  collision  broadening  effect  enhances  consistently  the  impact 
ionization,  although  the  differences  are  almost  within  the  error  of  the 
simulation.  Figure  6.7  shows  the  analogous  results  for  Si.  Again,  the 
collision  broadening  effect  enhances  the  impact  ionization.  This  can  be 
understood  because  the  collision  broadened  spectral  density  allows  some 
finite  probability  for  an  electron  in  a  momentum  state  k  to  have  an  energy 
greater  than  E(k).  Consider  the  following  "thought"  experiment: 

Suppose  we  put  an  electron  in  a  state  E(k)  in  a  solid  where  inter¬ 
action  between  the  electron  and  the  lattice  can  be  turned  off.  When 
there  is  no  interaction,  the  electron  stays  in  the  same  state  as  time 
evolves.  We  then  open  "window"  such  that  electrons  with  energy  E(k) 
escape  out  of  the  solid  with  the  following  probability: 

fl  for  E  (k)  >  E 

P  =  \  ° 

esc  I 

for  E(k)  <  E 

o 


0 


(6.44) 
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Figure  6.5:  The  imaginary  part  of  the  self-energy  for 
different  coupling  constants. 
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Now  we  turn  the  interaction  on  and  the  electron  starts  to  behave  like 
a  quasi  particle  with  a  broadened  energy  spectrum  as  time  evolves.  If 
we  assume  the  same  collision  broadening  7  for  all  the  states  and  neglect 
the  shifts  in  energy  caused  by  the  real  part  of  the  self-energy,  the 
escaping  probability  for  an  electron  in  a  state  E(k)  becomes  a  folded 
integral 
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esc 
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esc 
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(6.45) 


where  A(E' ,E) ,  the  normalized  spectral  density  function  from  Equation 
6.25,  is  expressed  as 
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Equation  6.44  can  be  analytically  integrated  to  give 
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The  collision  broadened  probability  is  plotted  in  Figure  6.8.  We  see 
that  the  probability  is  smeared-out  around  E^,  then  the  effect  of  the 
collision  broadening  is  to  smear  the  sharp  edge  of  the  barrier.  In  a 
sense,  the  barrier  height  is  effectively  reduced  for  the  electrons  re¬ 
siding  in  the  energy  range  below  the  barrier  height.  This  justifies  the 
argument  in  Chapter  5  that  the  emission  barrier  at  the  Si-SiO.,  interface 
should  be  further  reduced  due  to  the  collision  broadening  effect. 


I 
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The  "experiment"  described  above  is  by  no  means  rigorous.  It  des¬ 
cribes,  at  most,  the  qualitative  nature  of  the  collision  broadening.  A 
better  understanding  of  the  whole  problem  would  require  a  full-blown 
many  body  quantum  mechanical  treatment  which  should  be  a  full  research 
project  in  its  own  right. 

Finally,  we  show  in  Figure  6.9  the  intra-collisional  field  effect 
(ICFE)  for  the  steady  state  high  field  transport  in  Si  and  GaAs.  ICFE 
is  defined  here  to  be  the  average  ratio  of  the  energy  gained  in  the 
collision  duration  to  the  collision  broadening  half  width  as  given  in 
Equation  6.17.  We  see  that  the  percentage  effect  is  very  small.  As  dis¬ 
cussed  in  Section  6.3.2(ii),  it  is  simply  because  the  distribution  is  so 
heated  up  by  the  high  electric  field  that  the  electrons  are  mostly  re¬ 
siding  in  the  high  energy  range  where  the  ICFE  is  negligible.  This  ef¬ 
fect  might  be  important  in  the  initial  high  field  transient  regime  as 
discussed  by  Barker  [86].  Notice,  however,  that  the  formula  for  the  ICFE 
was  derived  by  an  effective  mass  approximation.  The  validity  of  the 
formula  in  high  energies  needs  further  justification. 

6.5  Summary 

The  quantum  transport  theory  has  been  in  existence  for  many  years 
[78-86].  Despite  the  many  efforts  of  trying  to  construct  a  quantum  trans¬ 
port  equation,  which  describes  transport  physics  better,  they  always  end 
up  facing  an  insolvable  transport  equation.  This  seriously  limits  the 
progress  in  the  quantitative  understanding  of  the  transport  problem.  In 
this  work,  we  have  discussed  the  assumptions  of  the  classical  Boltzmann 
equation  to  the  quantum  self-energy  effect  and  proposed  a  quantum  Monte 
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Carlo  method.  Our  quantum  Monte  Carlo  method  differs  from  the  classical 
one  in  that  the  quasi-particle  charac ter ist ics  have  been  taken  care  of 
in  the  simulation.  A  rigorous  justification  of  the  method  proposed  will 
require  further  investigation.  It  is  hoped  that  the  quantum  Monte  Carlo 
method  developed  in  this  work  will  serve  as  a  starting  basis  for  future 
research. 

As  for  the  transport  parameters,  the  deformation  potentials  which 
crucially  determine  the  scattering  rates  have  been  treated  essentially 
as  adjustable  parameters.  A  first  principal  theory  for  calculating  the 
deformation  potential  is  of  utmost  importance. 

It  has  been  demonstrated  in  this  thesis  work  that  the  Monte  Carlo 
simulation  can  be  used  to  treat  a  large  variety  of  transport  problems  in¬ 
accessible  by  any  other  means.  Future  challenges  in  improving  the  Monte 
Carlo  method  will  be  to  include  the  electron-electron  scattering;  to 
refine  the  band  structure  interpolation  scheme;  to  calculate  the  scattering 
rates,  including  the  quantum  effects;  etc.  The  Monte  Carlo  method  is 
powerful,  flexible,  and  not  limited  by  itself.  The  only  limitations  to 
this  method  are  the  capacity  and  speed  of  the  computer  and  the  lack  of  the 
fundamental  understanding  of  the  basic  transport  parameters. 
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APPENDIX  1 

DETERMINATION  OF  THE  CARRIER  FREE  FLIGHT  TIME 
IN  A  MONTE  CARLO  SIMULATION 

In  a  Monte  Carlo  simulation,  the  determination  of  the  free  flight 
time  is  one  of  the  most  important  objectives.  A  better  method  improves 
the  speed  of  the  program  and  the  accuracy  of  the  simulation.  The  free 
flight  time  is  the  time  that  the  electron  drifts  freely  between  succes¬ 
sive  collisions.  The  lengths  of  the  free  flight  times  follow  a  distri¬ 
bution  directly  related  to  the  scattering  rates  in  the  semiconductor. 

In  this  appendix,  we  derive  this  distribution  and  discuss  how  this  dis¬ 
tribution  can  be  properly  implemented  in  the  Monte  Carlo  simulation. 

Suppose  that  P(t^)  is  the  probability  that  an  electron  drifts 
freely  from  time  t  =  0  to  time  t  =  t ^ .  The  probability  for  the  electron 
not  to  be  scattered  in  the  time  interval  t^<t<t^+At  is  proportional  to 
the  total  scattering  rate  of  the  electron  and  the  length  of  the  interval 
as  At/x(t^).  So  the  probability  of  the  electron  not  getting  scattered 
in  this  interval  is  simply  l-At/x(t^).  It  follows  immediately  that  the 
probability  of  the  electron  drifting  freely  from  time  t=0  to  time  t=t^+At 
can  be  written  as 

P(tL+At)  -  P(t1)(l  -  T777)  •  (  Al.l) 

If  At  is  an  infinitesimal  increment  of  time,  which  then  makes  Equation 
Al.l  exact,  P(t)  is  solved  as 
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P(t  )  =  exp 


(A1.2) 


where  P(0)  is  assumed  to  be  one.  If  T  is  a  constant,  then  P(t)  is 
simply  a  Poisson  distribution  and  the  mean  of  which  is  just  T,  a  con¬ 
stant  scattering  mean  free  time.  Unfortunately,  T  is  a  complicated 
function  of  the  crystal  momenta  and  is  usually  simplified  to  be  a 
function  of  the  energy  which  implies  the  assumption  of  isotropic  scatter¬ 
ing. 

What  is  required  in  a  Monte  Carlo  simulation  is  a  means  to  deter¬ 
mine  the  time  interval  T  in  which  the  electron  drifts  freely,  and  is 
scattered  at  the  end  of  the  interval.  From  Equation  A1.2,  and  knowing 
that  the  probability  per  unit  time  of  having  one  scattering  event  at 
t=t-^  is  just  the  scattering  rate  at  t=t^ ,  the  distribution  for  T  can  be 
expressed  as 


P  (T)  =  P(t) 
sc 
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x(T) 


Ten  exp 


T 

r  dt ' 
T(t') 
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(Al. 3) 


By  generating  a  random  number  r  uniformly  distributed  between  0  and  1, 

T  can  be  obtained  by  solving 
T 

r  =  P  (t)  dt 
J  sc 
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that  the  probability  per  unit  time  of  having  one  scattering  event  at 
t=t^  is  just  the  scattering  rate  at  t=t^,  the  distribution  for  T  can  be 
expressed  as 


By  generating  a  random  number  r  uniformly  distributed  between  0  and  1, 

T  can  be  obtained  by  solving 
T 

r  =  P  (t)  dt 
I  sc 


0 
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0 


exp 


dc' 

.  T(t’) 

0 


dt 


(A1.4) 


Since  there  is  no  difference  in  probability  for  r  or  (1-r)  to  occur. 
Equation  (Ai.4)  is  sometimes  written  as 


r 


exp 


I 

|  r  dt 

'  J  T(t) 

L  0  _ 


(A1.5) 


As  mentioned  earlier,  t(E(k(t)))  is  a  complicated  function,  and  its 
dependence  on  time  is  through  the  wave  vector,  the  equation  of  motion, 
and  hence  the  band  structure.  There  is  no  simple  way  of  solving  Equation 
A1 . 5 .  One  has  to  invent  some  mathematical  device  to  solve  the  equation 
or  to  find  a  way  to  get  around  it.  We  describe  two  methods  which  have 
been  commonly  used  to  solve  this  problem. 

As  discussed  in  Reference  16,  one  can  adopt  a  pseudo  scattering 
mechanism  to  simplify  Equation  A1.5.  This  is  commonly  called  a  self¬ 
scattering  method.  It  is  a  mathematical  device  without  physical  sig¬ 
nificance  and  the  idea  of  which  is  summarized  as  follows. 


One  includes  in  Che  simulation  an  extra  scattering  mechanism  of 


the  form 


S  (k,k' )  = 
ps 


(E (k) ) 


•  (k  -  k')  , 


(A1.6) 


where  S  (k,k' )  is  the  transition  probability  from  state  k  to  k' ,  and 
ps 

— ► 

1/x  (E ( k) )  is  the  scattering  rate  for  the  mechanism.  Because  the 
ps 

6-function  conserves  the  wave  vector  of  the  carrier,  i.e.,  the  carrier 
is  exchanged  by  itself,  one  sees  that  this  mechanism  is  physically  in¬ 
significant.  If  the  scattering  strength  is  chosen  as 


(E)  ’ 


(Al.  71 


where  1/x^,  is  a  constant  and  is  now  the  total  scattering  rate,  one  sees 
that  Equation  A1.5  can  be  easily  solved  as 


r  =  exp 


J  T  (t)  T(t) 
0  ^  ps 


TTo] 


=  exp 


!  —  dt 


exp  - 


(  A 1 . 8  > 


Thus,  the  free  flight  time  is  simply  expressed  as 
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T  =  TT  ln  r  '  (A1 . 9) 

In  an  actual  simulation,  when  self-scattering  occurs,  one  simply 
lets  the  carrier  continue  in  its  present  momentum  state.  The  only  draw¬ 
back  with  the  self-scattering  scheme  is  that  the  total  constant  scattering 
rate  has  to  be  greater  than  the  maximum  true  scattering  rate  in  order  to 
make  the  self-scattering  rate  always  positive.  Thus,  the  self-scattering 
rate,  which  is  complementary  to  the  true  total  scattering  rate,  is  ex¬ 
tremely  high  at  low  energies.  This  makes  the  simulation  inefficient  be¬ 
cause  most  of  the  CPU  time  is  spent  on  counting  self-scattering  events 
which  have  no  physical  significance.  In  order  to  improve  the  efficiency  of 
the  simulation,  one  can  divide  the  energy  space  into  regions  and  assign 
different  values  of  x  to  different  regions  to  optimize  the  simulation. 

But  when  carriers  propagate  from  one  region  to  the  other,  one  again  deals 
with  complicated  mathematics  which  is  simply  not  practical  if  a  realistic 
band  structure  is  considered.  If  one  considers  a  simple  parabolic  band, 
this  may  be  a  way  to  go.  The  details  for  the  discussion  can  be  found  in 
[17]  and  are  omitted  here. 

The  second  method  described  below  was  what  has  been  used  in  our  Monte 
Carlo  simulation.  The  central  idea  is  to  bypass  completely  the  complicated 
integral  form  of  Equation  A1.4  and  go  back  to  its  differential  form. 

As  discussed  at  the  beginning  of  the  appendix,  the  probability  of 
a  carrier  being  scattered  in  the  short  interval  t^'t^t^+it  can  be  written 
as 


P 

sc 


(t. 


;At) 


)  ‘ 


(Al. 10) 
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This  holds  true  as  long  as  at  is  short  enough  that  i(t)  remains  almost 
a  constant  in  the  interval.  In  the  simulation,  if  we  let  the  carrier 
drift  during  a  small  increment  of  time  such  that  Equation  A1.10  holds 
true,  then  random  numbers  can  be  generated  to  determine  whether  the 
carrier  is  to  be  scattered  in  the  next  time  increment. 

The  actual  Monte  Carlo  procedure  is  described  as  follows.  We  first 
find  the  total  scattering  rate  correspond ing  to  the  present  state  of  the 
carrier.  We  then  let  it  drift  freely  for  one  tenth  of  the  mean  free 
scattering  time,  which  is  the  reciprocal  of  the  total  scattering  rate.  A 
random  number  r,  uniformiv  distributed  between  0  and  I,  is  generated  to 
represent  the  probability  or  scattering  occurring  in  the  next  time  incre¬ 
ment  At.  If  t!ie  generated  random  number  r  satisfies  the  relation 


the  carrier  is  scattered.  This  process  is  conceptually  simple  and  is 
exact  if  At  is  an  infinitesimal  increment.  In  a  practical  simulation,  one 
can  not  afford  to  use  too  small  a  time  increment  because  it  makes  the  pro¬ 
gram  extremely  slow.  As  long  as  At  is  chosen  to  be  about  one  tenth  the  mean 
free  scattering  time,  one  finds  that  both  the  speed  of  the  program  and  the 
accuracy  of  the  calculation  are  acceptable.  So,  in  our  simulation,  At  is 
always  chosen  to  be  about  one  tenth  of  the  mean  free  scattering  time  cor¬ 
responding  to  the  present  state  of  the  carrier. 
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APPENDIX  2 

SUMMARY  OF  SCATTERING  MECHANISMS  IN  SILICON 

In  this  appendix,  we  summarize  the  different  scattering  mechanisms 
that  have  been  included  in  our  Monte  Carlo  model  for  silicon.  As  dis¬ 
cussed  in  Section  2.2,  the  low  energy  scattering  rates  are  calculated  by 
the  conventional  method:  the  Golden  rule,  the  Born  approximation  and  the 
effective  mass  density-of-states. 

The  analytic  band  structure  for  low  energy  near  the  band  edge  of  an 
•*> 

electron  in  the  state  k  in  the  ith  valley  is  given  by 


vr(i),tu  _ 

Y  (c.  (k)  )  =  — 


,?  £(i).  2  t  -*(i)  2  | 

u-k0  )  L  U-k0  )  t| 


<  A  2  .  1 ) 


1 


and 


V(E)  =  E(l-KtE) 


(A- . 2 ) 


-►  ( i ) 

where  :t  is  the  nonparabol ic ity  factor;  k^  is  the  wave  vector  ot  the 
minimum  of  the  ith  valley;  the  subscripts  1  and  t.  stand  for  longitudinal 
and  transverse  components  with  respect  to  the  symmetry  axis  of  the  valley; 
m^  and  m  are  the  longitudinal  and  transverse  effective  masses  of  the 
electrons.  The  density-of-states  effective  mass  corresponding  to  the 
elliptic  energy  surface  is  given  by 


( A2 .  3) 


14  2 


m 


D 


(mimt 


1/3 


In  all  calculations,  the  overlap  integrals  are  assumed  to  be  one  and 
the  squared  matrix  elements  are  always  expressed  in  terms  of  the  assumed 
constant  deformation  potentials.  The  values  of  the  deformation  potential 
constants  are  determined  mainly  through  fits  to  experimental  transport 
data,  for  example,  drift  velocities,  impact  ionization  coefficients,  etc. 
We  follow  closely  the  work  of  Canali  et  al.  [18]  and  the  material  param¬ 
eters  are  listed  in  Appendix  3. 


2 . 1  Intravallev  acoustic  scattering 

As  discussed  in  [18],  the  angular  dependence  of  the  matrix  element 
for  acoustic  phonon  scattering  can  be  averaged  into  the  squared  matrix 
element  given  by  [24] 


H' 


9 


E“tfq 

2Vpv 

s 


In  +  i 
l  <i 


(  Al .  4) 


where  E,  is  the  deformation-potential  constant,  V  and  p  are 

ind  the  densitv  of  the  semiconductor,  v  is  the  velocity  of 

s 

N  is  the  phonon  distribution;  N  or  N  +1  must  be  taken  for 

q  q  q 

emission,  respectively.  By  application  of  the  Golden  rule, 
-lament  in  Equation  A2 . 4  yields  the  scattering  rate 


the  volume 
sound,  and 
absorpt ion 
the  matrix 


or 


p 
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P  (q)d  q 
ac 


Eb 


i 

8ir~pv 


[  *  6  (  E  ( k+q  )  - 


E(k)  +  fiqv  )d  q  , 


(A2. 5) 


N  +  1 


where  the  upper  or  lower  sign  must  be  taken  for  absorption  or  emission, 
respectively.  In  order  to  simplify  the  mathematics,  the  nonparabolicity 
factor  a  is  put  to  zero  for  acoustic  scattering.  So  by  performing  the 
Herring  and  Vogt  [104]  transformation  with  oc=0.  Equation  A2.1  becomes 


2  *  ’>  ( i ) 

.  (i)  _  f  k 


(A2.6) 


-*(i)  -  -  (i) 

where  k  is  the  transformed  wave  vector  of  k-k^  in  the  ith  valley, 

and  is  the  free  mass  of  the  electrons.  Transforming  to  the  starred 

space,  one  again  approximates 


r 

1/2 

(  \ 

★ 

ml  2  *  ^  mt  .  2„*i 

^  /’< 

q  =  q 

1 

! —  cos  u  4*  —  sin  v  | 

imo  mo  J 

=  q 

— 

1/2 


(A2.7) 


Vc 

where  0  is  the  angle  between  q  and  the  principal  axis  of  the  valley, 
and  trip  the  density  of  states  effective  mass.  The  argument  of  the 
energy  conservation  6-function  then  becomes 


„2  *2  „2  ,  ,  .  fm 

tfq  ,  -R  ,  *  *  -  ^  *  i  D 

— r-1 —  +  —  k  q  cosy  +  fiq  I 

2m0  ~  mQ  ^Qj 


1/2 


(A2 . 6 ) 
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where  y  is  the  angle  between  q  and  k  .  As  in  the  standard  procedure 
[24],  Equation  A2.7  is  put  to  zero  and  the  condition  |cosy|  <  1  yields 


the  limits  for  q  .  It  is  more  convenient  to  express  these  limits  by  a 
dimensionless  variable 


•fiqv 

s 

f  *  > 
•fiq  v 

s 

kT 

kT 

0J 

71/2 


( A2 . 9) 


These  limits  are  given  in  Table  A2.1  for  different  cases,  where 
*  1  2 

E  =  —  n>QVs*  The  transition  rate  is  given  by 


P 

ac 


(x)a^x 


Nq(x) 


N  (x)  +  ll 

l  9  J 


2 

x“dx 


f(x)dx  , 


(A2.10) 


where 


EimD1 

->  , 

1  1 

fH  1 

L  U*2  . 
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^Es 

(A2.ll) 


At  this  point,  one  must  introduce  an  expression  for  N  (x).  A 

q 

suitable  approximation  is  given  by  the  following  truncated  Laurent 
expansion  [ 18] : 


1 

720 


3.5 


N  (x) 

q 


0 


X  >  3.5 


(A2.12) 
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Table  A2.1  [18]  : 

Integration  limits  for  Equation  A2.9 

with  x  =  hqV  /kT. 
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E  <  E 


Absorption  < 


\  a2 


4E 


4E 


*1/2 
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No  emission 


Absorption  < 


x,  =  0 


4E 


*1/2 

t^<e1/2-;i/2> 


Emission 


t  \  -  o 


4E 


*1/2 


kT 


(E1/2  -  E  *1/2  ) 
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The  integration  of  the  scattering  rate  in  Equation  A2.9  can  then 
be  performed  very  easily  as 


f (x)dx  ,  (A2.13) 


ac 


where  x^  and  x0  are  listed  in  Table  A2.1.  The  resulting  scattering  rates 
are  given  in  Table  A2.2  and  plotted  in  Figure  A2.1. 

*  _5 

For  room  temperature,  E^  =  7.46x10  eV  and  x  corresponds  to  an 

—  k 

energy  of  about  7  eV.  The  cases  for  x^  >  x  and  E  <  E  are  never  really 
used  in  actual  calculations.  For  lower  temperatures,  all  cases  have  to 
be  considered. 

The  energy  exchanged  in  the  acoustic  scattering  depends  on  the 

momentum  exchanged  during  the  process.  To  keep  track  of  both  energy  and 

momentum  exchanged  in  the  process  requires  an  analytic  form  of  the  band 

structure.  To  simplify  the  situation,  we  consider  an  average  phonon 

energy  as  the  energy  exchanged  during  the  acoustic  scattering  processes. 

The  average  phonon  energy,  which  is  a  function  of  the  electron  energy,  is 

obtained  from  the  distribution  f(x)  of  the  phonon  momentum  x  =  -ffqv  ,'kT 

s 

given  in  Equation  A2.9.  For  a  given  electron  energy  x?(E)  (see  Table 
A2.1),  the  corresponding  average  phonon  energy  can  be  obtained  as 


xf  (x )dx 


f (x)dx 


<x> 


(A2. 14) 


Table  A2.2  (18]: 
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or 


<hoj  >  =  kT<x>  .  (A2.15) 

q 

We  plot  in  Figure  A2.2  the  average  phonon  energies  for  both  emission 
and  absorption  at  room  temperature.  Notice  that  at  high  electron 
energies,  the  energies  exchanged  in  the  scattering  process  are  not  negli¬ 
gible  which  makes  the  usual  assumption  of  equipartition  questionable. 

2 . 2  Intervalley  scattering 

The  scattering  mechanism  is  treated  here  in  the  traditional  way 
[24].  For  equivalent  intervalley  scattering,  the  total  scattering  rate 
is  given  by 
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where 
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Figure  A2.2:  The  average  acoustic  phonon  energy  as  a  function 
of  the  electron  energy  for  both  emission  and 
absorpt ion . 
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and  uk  are  the  coupling  constant  and  the  phonon  energy  (assumed  con¬ 
stant)  of  the  considered  ith  mechanism;  is  the  number  of  possible 
final  valleys.  For  nonequivalent  intervalley  scattering,  it  is  essentially 
given  by  the  same  formula  except  that  E'  is  referred  to  the  proper  valley 
minima.  Details  of  the  formula  can  be  found  in  appendix  of  Reference  25 
and  is  hence  omitted  here. 

The  scattering  rates  for  f  and  g  scatterings  are  plotted  in  Figure 
A2.3  respectively. 
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Figure  A2.3:  Intervalley  scattering  rates  in  silicon. 
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APPENDIX  3 

MATERIAL  PARAMETERS  FOR  SILICON 

Bulk  Material  Parameters 


Lattice  constant 

5.43 

O 

A 

Density 

2.  329 

g/ cm' 

Dielectric  constant 

11.7 

Sound  velocity 

9.04xl05 

cm/  s 

X-Valley 

Effective  masses 

transverse 

0.19 

mo 

longitudinal 

0.916 

rao 

Nonparabolicitv 

0.5 

eV': 

Acoustic  deformation 
po  tent ial 

9.5 

eV 

L-Vallev 

Effective  masses 

transverse 

0.12 

m0 

longitudinal 

1.59 

mQ 
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X-X  Intervaliey  Scattering 


Phonon 

temperature (K) 

Deformation 
potent ial (eV/cm) 

Scattering 

220 

3xl07 

f 

550 

2xl08 

f 

685 

2xl08 

f 

140 

5xl07 

g 

215 

8xl07 

g 

720 

9 

1. lxlO7 

g 

X-L  Intervalley  Scattering 


Phonon 

temperature (K) 

672 

634 

480 


Deformation 
potent ial (eV/cm) 

2xl08 

2xl08 

2xl08 

2xl08 


197 
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